diff options
Diffstat (limited to 'src/backend/utils/adt/geo_ops.c')
| -rw-r--r-- | src/backend/utils/adt/geo_ops.c | 5569 |
1 files changed, 5569 insertions, 0 deletions
diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c new file mode 100644 index 0000000..9484dbc --- /dev/null +++ b/src/backend/utils/adt/geo_ops.c @@ -0,0 +1,5569 @@ +/*------------------------------------------------------------------------- + * + * geo_ops.c + * 2D geometric operations + * + * This module implements the geometric functions and operators. The + * geometric types are (from simple to more complicated): + * + * - point + * - line + * - line segment + * - box + * - circle + * - polygon + * + * Portions Copyright (c) 1996-2021, PostgreSQL Global Development Group + * Portions Copyright (c) 1994, Regents of the University of California + * + * + * IDENTIFICATION + * src/backend/utils/adt/geo_ops.c + * + *------------------------------------------------------------------------- + */ +#include "postgres.h" + +#include <math.h> +#include <limits.h> +#include <float.h> +#include <ctype.h> + +#include "libpq/pqformat.h" +#include "miscadmin.h" +#include "utils/float.h" +#include "utils/fmgrprotos.h" +#include "utils/geo_decls.h" + +/* + * * Type constructors have this form: + * void type_construct(Type *result, ...); + * + * * Operators commonly have signatures such as + * void type1_operator_type2(Type *result, Type1 *obj1, Type2 *obj2); + * + * Common operators are: + * * Intersection point: + * bool type1_interpt_type2(Point *result, Type1 *obj1, Type2 *obj2); + * Return whether the two objects intersect. If *result is not NULL, + * it is set to the intersection point. + * + * * Containment: + * bool type1_contain_type2(Type1 *obj1, Type2 *obj2); + * Return whether obj1 contains obj2. + * bool type1_contain_type2(Type1 *contains_obj, Type1 *contained_obj); + * Return whether obj1 contains obj2 (used when types are the same) + * + * * Distance of closest point in or on obj1 to obj2: + * float8 type1_closept_type2(Point *result, Type1 *obj1, Type2 *obj2); + * Returns the shortest distance between two objects. If *result is not + * NULL, it is set to the closest point in or on obj1 to obj2. + * + * These functions may be used to implement multiple SQL-level operators. For + * example, determining whether two lines are parallel is done by checking + * whether they don't intersect. + */ + +/* + * Internal routines + */ + +enum path_delim +{ + PATH_NONE, PATH_OPEN, PATH_CLOSED +}; + +/* Routines for points */ +static inline void point_construct(Point *result, float8 x, float8 y); +static inline void point_add_point(Point *result, Point *pt1, Point *pt2); +static inline void point_sub_point(Point *result, Point *pt1, Point *pt2); +static inline void point_mul_point(Point *result, Point *pt1, Point *pt2); +static inline void point_div_point(Point *result, Point *pt1, Point *pt2); +static inline bool point_eq_point(Point *pt1, Point *pt2); +static inline float8 point_dt(Point *pt1, Point *pt2); +static inline float8 point_sl(Point *pt1, Point *pt2); +static int point_inside(Point *p, int npts, Point *plist); + +/* Routines for lines */ +static inline void line_construct(LINE *result, Point *pt, float8 m); +static inline float8 line_sl(LINE *line); +static inline float8 line_invsl(LINE *line); +static bool line_interpt_line(Point *result, LINE *l1, LINE *l2); +static bool line_contain_point(LINE *line, Point *point); +static float8 line_closept_point(Point *result, LINE *line, Point *pt); + +/* Routines for line segments */ +static inline void statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2); +static inline float8 lseg_sl(LSEG *lseg); +static inline float8 lseg_invsl(LSEG *lseg); +static bool lseg_interpt_line(Point *result, LSEG *lseg, LINE *line); +static bool lseg_interpt_lseg(Point *result, LSEG *l1, LSEG *l2); +static int lseg_crossing(float8 x, float8 y, float8 px, float8 py); +static bool lseg_contain_point(LSEG *lseg, Point *point); +static float8 lseg_closept_point(Point *result, LSEG *lseg, Point *pt); +static float8 lseg_closept_line(Point *result, LSEG *lseg, LINE *line); +static float8 lseg_closept_lseg(Point *result, LSEG *on_lseg, LSEG *to_lseg); + +/* Routines for boxes */ +static inline void box_construct(BOX *result, Point *pt1, Point *pt2); +static void box_cn(Point *center, BOX *box); +static bool box_ov(BOX *box1, BOX *box2); +static float8 box_ar(BOX *box); +static float8 box_ht(BOX *box); +static float8 box_wd(BOX *box); +static bool box_contain_point(BOX *box, Point *point); +static bool box_contain_box(BOX *contains_box, BOX *contained_box); +static bool box_contain_lseg(BOX *box, LSEG *lseg); +static bool box_interpt_lseg(Point *result, BOX *box, LSEG *lseg); +static float8 box_closept_point(Point *result, BOX *box, Point *point); +static float8 box_closept_lseg(Point *result, BOX *box, LSEG *lseg); + +/* Routines for circles */ +static float8 circle_ar(CIRCLE *circle); + +/* Routines for polygons */ +static void make_bound_box(POLYGON *poly); +static void poly_to_circle(CIRCLE *result, POLYGON *poly); +static bool lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start); +static bool poly_contain_poly(POLYGON *contains_poly, POLYGON *contained_poly); +static bool plist_same(int npts, Point *p1, Point *p2); +static float8 dist_ppoly_internal(Point *pt, POLYGON *poly); + +/* Routines for encoding and decoding */ +static float8 single_decode(char *num, char **endptr_p, + const char *type_name, const char *orig_string); +static void single_encode(float8 x, StringInfo str); +static void pair_decode(char *str, float8 *x, float8 *y, char **endptr_p, + const char *type_name, const char *orig_string); +static void pair_encode(float8 x, float8 y, StringInfo str); +static int pair_count(char *s, char delim); +static void path_decode(char *str, bool opentype, int npts, Point *p, + bool *isopen, char **endptr_p, + const char *type_name, const char *orig_string); +static char *path_encode(enum path_delim path_delim, int npts, Point *pt); + + +/* + * Delimiters for input and output strings. + * LDELIM, RDELIM, and DELIM are left, right, and separator delimiters, respectively. + * LDELIM_EP, RDELIM_EP are left and right delimiters for paths with endpoints. + */ + +#define LDELIM '(' +#define RDELIM ')' +#define DELIM ',' +#define LDELIM_EP '[' +#define RDELIM_EP ']' +#define LDELIM_C '<' +#define RDELIM_C '>' +#define LDELIM_L '{' +#define RDELIM_L '}' + + +/* + * Geometric data types are composed of points. + * This code tries to support a common format throughout the data types, + * to allow for more predictable usage and data type conversion. + * The fundamental unit is the point. Other units are line segments, + * open paths, boxes, closed paths, and polygons (which should be considered + * non-intersecting closed paths). + * + * Data representation is as follows: + * point: (x,y) + * line segment: [(x1,y1),(x2,y2)] + * box: (x1,y1),(x2,y2) + * open path: [(x1,y1),...,(xn,yn)] + * closed path: ((x1,y1),...,(xn,yn)) + * polygon: ((x1,y1),...,(xn,yn)) + * + * For boxes, the points are opposite corners with the first point at the top right. + * For closed paths and polygons, the points should be reordered to allow + * fast and correct equality comparisons. + * + * XXX perhaps points in complex shapes should be reordered internally + * to allow faster internal operations, but should keep track of input order + * and restore that order for text output - tgl 97/01/16 + */ + +static float8 +single_decode(char *num, char **endptr_p, + const char *type_name, const char *orig_string) +{ + return float8in_internal(num, endptr_p, type_name, orig_string); +} /* single_decode() */ + +static void +single_encode(float8 x, StringInfo str) +{ + char *xstr = float8out_internal(x); + + appendStringInfoString(str, xstr); + pfree(xstr); +} /* single_encode() */ + +static void +pair_decode(char *str, float8 *x, float8 *y, char **endptr_p, + const char *type_name, const char *orig_string) +{ + bool has_delim; + + while (isspace((unsigned char) *str)) + str++; + if ((has_delim = (*str == LDELIM))) + str++; + + *x = float8in_internal(str, &str, type_name, orig_string); + + if (*str++ != DELIM) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + type_name, orig_string))); + + *y = float8in_internal(str, &str, type_name, orig_string); + + if (has_delim) + { + if (*str++ != RDELIM) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + type_name, orig_string))); + while (isspace((unsigned char) *str)) + str++; + } + + /* report stopping point if wanted, else complain if not end of string */ + if (endptr_p) + *endptr_p = str; + else if (*str != '\0') + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + type_name, orig_string))); +} + +static void +pair_encode(float8 x, float8 y, StringInfo str) +{ + char *xstr = float8out_internal(x); + char *ystr = float8out_internal(y); + + appendStringInfo(str, "%s,%s", xstr, ystr); + pfree(xstr); + pfree(ystr); +} + +static void +path_decode(char *str, bool opentype, int npts, Point *p, + bool *isopen, char **endptr_p, + const char *type_name, const char *orig_string) +{ + int depth = 0; + char *cp; + int i; + + while (isspace((unsigned char) *str)) + str++; + if ((*isopen = (*str == LDELIM_EP))) + { + /* no open delimiter allowed? */ + if (!opentype) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + type_name, orig_string))); + depth++; + str++; + } + else if (*str == LDELIM) + { + cp = (str + 1); + while (isspace((unsigned char) *cp)) + cp++; + if (*cp == LDELIM) + { + depth++; + str = cp; + } + else if (strrchr(str, LDELIM) == str) + { + depth++; + str = cp; + } + } + + for (i = 0; i < npts; i++) + { + pair_decode(str, &(p->x), &(p->y), &str, type_name, orig_string); + if (*str == DELIM) + str++; + p++; + } + + while (depth > 0) + { + if (*str == RDELIM || (*str == RDELIM_EP && *isopen && depth == 1)) + { + depth--; + str++; + while (isspace((unsigned char) *str)) + str++; + } + else + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + type_name, orig_string))); + } + + /* report stopping point if wanted, else complain if not end of string */ + if (endptr_p) + *endptr_p = str; + else if (*str != '\0') + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + type_name, orig_string))); +} /* path_decode() */ + +static char * +path_encode(enum path_delim path_delim, int npts, Point *pt) +{ + StringInfoData str; + int i; + + initStringInfo(&str); + + switch (path_delim) + { + case PATH_CLOSED: + appendStringInfoChar(&str, LDELIM); + break; + case PATH_OPEN: + appendStringInfoChar(&str, LDELIM_EP); + break; + case PATH_NONE: + break; + } + + for (i = 0; i < npts; i++) + { + if (i > 0) + appendStringInfoChar(&str, DELIM); + appendStringInfoChar(&str, LDELIM); + pair_encode(pt->x, pt->y, &str); + appendStringInfoChar(&str, RDELIM); + pt++; + } + + switch (path_delim) + { + case PATH_CLOSED: + appendStringInfoChar(&str, RDELIM); + break; + case PATH_OPEN: + appendStringInfoChar(&str, RDELIM_EP); + break; + case PATH_NONE: + break; + } + + return str.data; +} /* path_encode() */ + +/*------------------------------------------------------------- + * pair_count - count the number of points + * allow the following notation: + * '((1,2),(3,4))' + * '(1,3,2,4)' + * require an odd number of delim characters in the string + *-------------------------------------------------------------*/ +static int +pair_count(char *s, char delim) +{ + int ndelim = 0; + + while ((s = strchr(s, delim)) != NULL) + { + ndelim++; + s++; + } + return (ndelim % 2) ? ((ndelim + 1) / 2) : -1; +} + + +/*********************************************************************** + ** + ** Routines for two-dimensional boxes. + ** + ***********************************************************************/ + +/*---------------------------------------------------------- + * Formatting and conversion routines. + *---------------------------------------------------------*/ + +/* box_in - convert a string to internal form. + * + * External format: (two corners of box) + * "(f8, f8), (f8, f8)" + * also supports the older style "(f8, f8, f8, f8)" + */ +Datum +box_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + BOX *box = (BOX *) palloc(sizeof(BOX)); + bool isopen; + float8 x, + y; + + path_decode(str, false, 2, &(box->high), &isopen, NULL, "box", str); + + /* reorder corners if necessary... */ + if (float8_lt(box->high.x, box->low.x)) + { + x = box->high.x; + box->high.x = box->low.x; + box->low.x = x; + } + if (float8_lt(box->high.y, box->low.y)) + { + y = box->high.y; + box->high.y = box->low.y; + box->low.y = y; + } + + PG_RETURN_BOX_P(box); +} + +/* box_out - convert a box to external form. + */ +Datum +box_out(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + + PG_RETURN_CSTRING(path_encode(PATH_NONE, 2, &(box->high))); +} + +/* + * box_recv - converts external binary format to box + */ +Datum +box_recv(PG_FUNCTION_ARGS) +{ + StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); + BOX *box; + float8 x, + y; + + box = (BOX *) palloc(sizeof(BOX)); + + box->high.x = pq_getmsgfloat8(buf); + box->high.y = pq_getmsgfloat8(buf); + box->low.x = pq_getmsgfloat8(buf); + box->low.y = pq_getmsgfloat8(buf); + + /* reorder corners if necessary... */ + if (float8_lt(box->high.x, box->low.x)) + { + x = box->high.x; + box->high.x = box->low.x; + box->low.x = x; + } + if (float8_lt(box->high.y, box->low.y)) + { + y = box->high.y; + box->high.y = box->low.y; + box->low.y = y; + } + + PG_RETURN_BOX_P(box); +} + +/* + * box_send - converts box to binary format + */ +Datum +box_send(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + StringInfoData buf; + + pq_begintypsend(&buf); + pq_sendfloat8(&buf, box->high.x); + pq_sendfloat8(&buf, box->high.y); + pq_sendfloat8(&buf, box->low.x); + pq_sendfloat8(&buf, box->low.y); + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); +} + + +/* box_construct - fill in a new box. + */ +static inline void +box_construct(BOX *result, Point *pt1, Point *pt2) +{ + if (float8_gt(pt1->x, pt2->x)) + { + result->high.x = pt1->x; + result->low.x = pt2->x; + } + else + { + result->high.x = pt2->x; + result->low.x = pt1->x; + } + if (float8_gt(pt1->y, pt2->y)) + { + result->high.y = pt1->y; + result->low.y = pt2->y; + } + else + { + result->high.y = pt2->y; + result->low.y = pt1->y; + } +} + + +/*---------------------------------------------------------- + * Relational operators for BOXes. + * <, >, <=, >=, and == are based on box area. + *---------------------------------------------------------*/ + +/* box_same - are two boxes identical? + */ +Datum +box_same(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(point_eq_point(&box1->high, &box2->high) && + point_eq_point(&box1->low, &box2->low)); +} + +/* box_overlap - does box1 overlap box2? + */ +Datum +box_overlap(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(box_ov(box1, box2)); +} + +static bool +box_ov(BOX *box1, BOX *box2) +{ + return (FPle(box1->low.x, box2->high.x) && + FPle(box2->low.x, box1->high.x) && + FPle(box1->low.y, box2->high.y) && + FPle(box2->low.y, box1->high.y)); +} + +/* box_left - is box1 strictly left of box2? + */ +Datum +box_left(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPlt(box1->high.x, box2->low.x)); +} + +/* box_overleft - is the right edge of box1 at or left of + * the right edge of box2? + * + * This is "less than or equal" for the end of a time range, + * when time ranges are stored as rectangles. + */ +Datum +box_overleft(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPle(box1->high.x, box2->high.x)); +} + +/* box_right - is box1 strictly right of box2? + */ +Datum +box_right(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPgt(box1->low.x, box2->high.x)); +} + +/* box_overright - is the left edge of box1 at or right of + * the left edge of box2? + * + * This is "greater than or equal" for time ranges, when time ranges + * are stored as rectangles. + */ +Datum +box_overright(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPge(box1->low.x, box2->low.x)); +} + +/* box_below - is box1 strictly below box2? + */ +Datum +box_below(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPlt(box1->high.y, box2->low.y)); +} + +/* box_overbelow - is the upper edge of box1 at or below + * the upper edge of box2? + */ +Datum +box_overbelow(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPle(box1->high.y, box2->high.y)); +} + +/* box_above - is box1 strictly above box2? + */ +Datum +box_above(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPgt(box1->low.y, box2->high.y)); +} + +/* box_overabove - is the lower edge of box1 at or above + * the lower edge of box2? + */ +Datum +box_overabove(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPge(box1->low.y, box2->low.y)); +} + +/* box_contained - is box1 contained by box2? + */ +Datum +box_contained(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(box_contain_box(box2, box1)); +} + +/* box_contain - does box1 contain box2? + */ +Datum +box_contain(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(box_contain_box(box1, box2)); +} + +/* + * Check whether the second box is in the first box or on its border + */ +static bool +box_contain_box(BOX *contains_box, BOX *contained_box) +{ + return FPge(contains_box->high.x, contained_box->high.x) && + FPle(contains_box->low.x, contained_box->low.x) && + FPge(contains_box->high.y, contained_box->high.y) && + FPle(contains_box->low.y, contained_box->low.y); +} + + +/* box_positionop - + * is box1 entirely {above,below} box2? + * + * box_below_eq and box_above_eq are obsolete versions that (probably + * erroneously) accept the equal-boundaries case. Since these are not + * in sync with the box_left and box_right code, they are deprecated and + * not supported in the PG 8.1 rtree operator class extension. + */ +Datum +box_below_eq(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPle(box1->high.y, box2->low.y)); +} + +Datum +box_above_eq(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPge(box1->low.y, box2->high.y)); +} + + +/* box_relop - is area(box1) relop area(box2), within + * our accuracy constraint? + */ +Datum +box_lt(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPlt(box_ar(box1), box_ar(box2))); +} + +Datum +box_gt(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPgt(box_ar(box1), box_ar(box2))); +} + +Datum +box_eq(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPeq(box_ar(box1), box_ar(box2))); +} + +Datum +box_le(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPle(box_ar(box1), box_ar(box2))); +} + +Datum +box_ge(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(FPge(box_ar(box1), box_ar(box2))); +} + + +/*---------------------------------------------------------- + * "Arithmetic" operators on boxes. + *---------------------------------------------------------*/ + +/* box_area - returns the area of the box. + */ +Datum +box_area(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + + PG_RETURN_FLOAT8(box_ar(box)); +} + + +/* box_width - returns the width of the box + * (horizontal magnitude). + */ +Datum +box_width(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + + PG_RETURN_FLOAT8(box_wd(box)); +} + + +/* box_height - returns the height of the box + * (vertical magnitude). + */ +Datum +box_height(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + + PG_RETURN_FLOAT8(box_ht(box)); +} + + +/* box_distance - returns the distance between the + * center points of two boxes. + */ +Datum +box_distance(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + Point a, + b; + + box_cn(&a, box1); + box_cn(&b, box2); + + PG_RETURN_FLOAT8(point_dt(&a, &b)); +} + + +/* box_center - returns the center point of the box. + */ +Datum +box_center(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + Point *result = (Point *) palloc(sizeof(Point)); + + box_cn(result, box); + + PG_RETURN_POINT_P(result); +} + + +/* box_ar - returns the area of the box. + */ +static float8 +box_ar(BOX *box) +{ + return float8_mul(box_wd(box), box_ht(box)); +} + + +/* box_cn - stores the centerpoint of the box into *center. + */ +static void +box_cn(Point *center, BOX *box) +{ + center->x = float8_div(float8_pl(box->high.x, box->low.x), 2.0); + center->y = float8_div(float8_pl(box->high.y, box->low.y), 2.0); +} + + +/* box_wd - returns the width (length) of the box + * (horizontal magnitude). + */ +static float8 +box_wd(BOX *box) +{ + return float8_mi(box->high.x, box->low.x); +} + + +/* box_ht - returns the height of the box + * (vertical magnitude). + */ +static float8 +box_ht(BOX *box) +{ + return float8_mi(box->high.y, box->low.y); +} + + +/*---------------------------------------------------------- + * Funky operations. + *---------------------------------------------------------*/ + +/* box_intersect - + * returns the overlapping portion of two boxes, + * or NULL if they do not intersect. + */ +Datum +box_intersect(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0); + BOX *box2 = PG_GETARG_BOX_P(1); + BOX *result; + + if (!box_ov(box1, box2)) + PG_RETURN_NULL(); + + result = (BOX *) palloc(sizeof(BOX)); + + result->high.x = float8_min(box1->high.x, box2->high.x); + result->low.x = float8_max(box1->low.x, box2->low.x); + result->high.y = float8_min(box1->high.y, box2->high.y); + result->low.y = float8_max(box1->low.y, box2->low.y); + + PG_RETURN_BOX_P(result); +} + + +/* box_diagonal - + * returns a line segment which happens to be the + * positive-slope diagonal of "box". + */ +Datum +box_diagonal(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + LSEG *result = (LSEG *) palloc(sizeof(LSEG)); + + statlseg_construct(result, &box->high, &box->low); + + PG_RETURN_LSEG_P(result); +} + +/*********************************************************************** + ** + ** Routines for 2D lines. + ** + ***********************************************************************/ + +static bool +line_decode(char *s, const char *str, LINE *line) +{ + /* s was already advanced over leading '{' */ + line->A = single_decode(s, &s, "line", str); + if (*s++ != DELIM) + return false; + line->B = single_decode(s, &s, "line", str); + if (*s++ != DELIM) + return false; + line->C = single_decode(s, &s, "line", str); + if (*s++ != RDELIM_L) + return false; + while (isspace((unsigned char) *s)) + s++; + if (*s != '\0') + return false; + return true; +} + +Datum +line_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + LINE *line = (LINE *) palloc(sizeof(LINE)); + LSEG lseg; + bool isopen; + char *s; + + s = str; + while (isspace((unsigned char) *s)) + s++; + if (*s == LDELIM_L) + { + if (!line_decode(s + 1, str, line)) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + "line", str))); + if (FPzero(line->A) && FPzero(line->B)) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid line specification: A and B cannot both be zero"))); + } + else + { + path_decode(s, true, 2, &lseg.p[0], &isopen, NULL, "line", str); + if (point_eq_point(&lseg.p[0], &lseg.p[1])) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid line specification: must be two distinct points"))); + line_construct(line, &lseg.p[0], lseg_sl(&lseg)); + } + + PG_RETURN_LINE_P(line); +} + + +Datum +line_out(PG_FUNCTION_ARGS) +{ + LINE *line = PG_GETARG_LINE_P(0); + char *astr = float8out_internal(line->A); + char *bstr = float8out_internal(line->B); + char *cstr = float8out_internal(line->C); + + PG_RETURN_CSTRING(psprintf("%c%s%c%s%c%s%c", LDELIM_L, astr, DELIM, bstr, + DELIM, cstr, RDELIM_L)); +} + +/* + * line_recv - converts external binary format to line + */ +Datum +line_recv(PG_FUNCTION_ARGS) +{ + StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); + LINE *line; + + line = (LINE *) palloc(sizeof(LINE)); + + line->A = pq_getmsgfloat8(buf); + line->B = pq_getmsgfloat8(buf); + line->C = pq_getmsgfloat8(buf); + + if (FPzero(line->A) && FPzero(line->B)) + ereport(ERROR, + (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION), + errmsg("invalid line specification: A and B cannot both be zero"))); + + PG_RETURN_LINE_P(line); +} + +/* + * line_send - converts line to binary format + */ +Datum +line_send(PG_FUNCTION_ARGS) +{ + LINE *line = PG_GETARG_LINE_P(0); + StringInfoData buf; + + pq_begintypsend(&buf); + pq_sendfloat8(&buf, line->A); + pq_sendfloat8(&buf, line->B); + pq_sendfloat8(&buf, line->C); + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); +} + + +/*---------------------------------------------------------- + * Conversion routines from one line formula to internal. + * Internal form: Ax+By+C=0 + *---------------------------------------------------------*/ + +/* + * Fill already-allocated LINE struct from the point and the slope + */ +static inline void +line_construct(LINE *result, Point *pt, float8 m) +{ + if (isinf(m)) + { + /* vertical - use "x = C" */ + result->A = -1.0; + result->B = 0.0; + result->C = pt->x; + } + else if (m == 0) + { + /* horizontal - use "y = C" */ + result->A = 0.0; + result->B = -1.0; + result->C = pt->y; + } + else + { + /* use "mx - y + yinter = 0" */ + result->A = m; + result->B = -1.0; + result->C = float8_mi(pt->y, float8_mul(m, pt->x)); + /* on some platforms, the preceding expression tends to produce -0 */ + if (result->C == 0.0) + result->C = 0.0; + } +} + +/* line_construct_pp() + * two points + */ +Datum +line_construct_pp(PG_FUNCTION_ARGS) +{ + Point *pt1 = PG_GETARG_POINT_P(0); + Point *pt2 = PG_GETARG_POINT_P(1); + LINE *result = (LINE *) palloc(sizeof(LINE)); + + if (point_eq_point(pt1, pt2)) + ereport(ERROR, + (errcode(ERRCODE_INVALID_PARAMETER_VALUE), + errmsg("invalid line specification: must be two distinct points"))); + + line_construct(result, pt1, point_sl(pt1, pt2)); + + PG_RETURN_LINE_P(result); +} + + +/*---------------------------------------------------------- + * Relative position routines. + *---------------------------------------------------------*/ + +Datum +line_intersect(PG_FUNCTION_ARGS) +{ + LINE *l1 = PG_GETARG_LINE_P(0); + LINE *l2 = PG_GETARG_LINE_P(1); + + PG_RETURN_BOOL(line_interpt_line(NULL, l1, l2)); +} + +Datum +line_parallel(PG_FUNCTION_ARGS) +{ + LINE *l1 = PG_GETARG_LINE_P(0); + LINE *l2 = PG_GETARG_LINE_P(1); + + PG_RETURN_BOOL(!line_interpt_line(NULL, l1, l2)); +} + +Datum +line_perp(PG_FUNCTION_ARGS) +{ + LINE *l1 = PG_GETARG_LINE_P(0); + LINE *l2 = PG_GETARG_LINE_P(1); + + if (FPzero(l1->A)) + PG_RETURN_BOOL(FPzero(l2->B)); + if (FPzero(l2->A)) + PG_RETURN_BOOL(FPzero(l1->B)); + if (FPzero(l1->B)) + PG_RETURN_BOOL(FPzero(l2->A)); + if (FPzero(l2->B)) + PG_RETURN_BOOL(FPzero(l1->A)); + + PG_RETURN_BOOL(FPeq(float8_div(float8_mul(l1->A, l2->A), + float8_mul(l1->B, l2->B)), -1.0)); +} + +Datum +line_vertical(PG_FUNCTION_ARGS) +{ + LINE *line = PG_GETARG_LINE_P(0); + + PG_RETURN_BOOL(FPzero(line->B)); +} + +Datum +line_horizontal(PG_FUNCTION_ARGS) +{ + LINE *line = PG_GETARG_LINE_P(0); + + PG_RETURN_BOOL(FPzero(line->A)); +} + + +/* + * Check whether the two lines are the same + */ +Datum +line_eq(PG_FUNCTION_ARGS) +{ + LINE *l1 = PG_GETARG_LINE_P(0); + LINE *l2 = PG_GETARG_LINE_P(1); + float8 ratio; + + /* If any NaNs are involved, insist on exact equality */ + if (unlikely(isnan(l1->A) || isnan(l1->B) || isnan(l1->C) || + isnan(l2->A) || isnan(l2->B) || isnan(l2->C))) + { + PG_RETURN_BOOL(float8_eq(l1->A, l2->A) && + float8_eq(l1->B, l2->B) && + float8_eq(l1->C, l2->C)); + } + + /* Otherwise, lines whose parameters are proportional are the same */ + if (!FPzero(l2->A)) + ratio = float8_div(l1->A, l2->A); + else if (!FPzero(l2->B)) + ratio = float8_div(l1->B, l2->B); + else if (!FPzero(l2->C)) + ratio = float8_div(l1->C, l2->C); + else + ratio = 1.0; + + PG_RETURN_BOOL(FPeq(l1->A, float8_mul(ratio, l2->A)) && + FPeq(l1->B, float8_mul(ratio, l2->B)) && + FPeq(l1->C, float8_mul(ratio, l2->C))); +} + + +/*---------------------------------------------------------- + * Line arithmetic routines. + *---------------------------------------------------------*/ + +/* + * Return slope of the line + */ +static inline float8 +line_sl(LINE *line) +{ + if (FPzero(line->A)) + return 0.0; + if (FPzero(line->B)) + return get_float8_infinity(); + return float8_div(line->A, -line->B); +} + + +/* + * Return inverse slope of the line + */ +static inline float8 +line_invsl(LINE *line) +{ + if (FPzero(line->A)) + return get_float8_infinity(); + if (FPzero(line->B)) + return 0.0; + return float8_div(line->B, line->A); +} + + +/* line_distance() + * Distance between two lines. + */ +Datum +line_distance(PG_FUNCTION_ARGS) +{ + LINE *l1 = PG_GETARG_LINE_P(0); + LINE *l2 = PG_GETARG_LINE_P(1); + float8 ratio; + + if (line_interpt_line(NULL, l1, l2)) /* intersecting? */ + PG_RETURN_FLOAT8(0.0); + + if (!FPzero(l1->A) && !isnan(l1->A) && !FPzero(l2->A) && !isnan(l2->A)) + ratio = float8_div(l1->A, l2->A); + else if (!FPzero(l1->B) && !isnan(l1->B) && !FPzero(l2->B) && !isnan(l2->B)) + ratio = float8_div(l1->B, l2->B); + else + ratio = 1.0; + + PG_RETURN_FLOAT8(float8_div(fabs(float8_mi(l1->C, + float8_mul(ratio, l2->C))), + HYPOT(l1->A, l1->B))); +} + +/* line_interpt() + * Point where two lines l1, l2 intersect (if any) + */ +Datum +line_interpt(PG_FUNCTION_ARGS) +{ + LINE *l1 = PG_GETARG_LINE_P(0); + LINE *l2 = PG_GETARG_LINE_P(1); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + if (!line_interpt_line(result, l1, l2)) + PG_RETURN_NULL(); + PG_RETURN_POINT_P(result); +} + +/* + * Internal version of line_interpt + * + * Return whether two lines intersect. If *result is not NULL, it is set to + * the intersection point. + * + * NOTE: If the lines are identical then we will find they are parallel + * and report "no intersection". This is a little weird, but since + * there's no *unique* intersection, maybe it's appropriate behavior. + * + * If the lines have NaN constants, we will return true, and the intersection + * point would have NaN coordinates. We shouldn't return false in this case + * because that would mean the lines are parallel. + */ +static bool +line_interpt_line(Point *result, LINE *l1, LINE *l2) +{ + float8 x, + y; + + if (!FPzero(l1->B)) + { + if (FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->B)))) + return false; + + x = float8_div(float8_mi(float8_mul(l1->B, l2->C), + float8_mul(l2->B, l1->C)), + float8_mi(float8_mul(l1->A, l2->B), + float8_mul(l2->A, l1->B))); + y = float8_div(-float8_pl(float8_mul(l1->A, x), l1->C), l1->B); + } + else if (!FPzero(l2->B)) + { + if (FPeq(l1->A, float8_mul(l2->A, float8_div(l1->B, l2->B)))) + return false; + + x = float8_div(float8_mi(float8_mul(l2->B, l1->C), + float8_mul(l1->B, l2->C)), + float8_mi(float8_mul(l2->A, l1->B), + float8_mul(l1->A, l2->B))); + y = float8_div(-float8_pl(float8_mul(l2->A, x), l2->C), l2->B); + } + else + return false; + + /* On some platforms, the preceding expressions tend to produce -0. */ + if (x == 0.0) + x = 0.0; + if (y == 0.0) + y = 0.0; + + if (result != NULL) + point_construct(result, x, y); + + return true; +} + + +/*********************************************************************** + ** + ** Routines for 2D paths (sequences of line segments, also + ** called `polylines'). + ** + ** This is not a general package for geometric paths, + ** which of course include polygons; the emphasis here + ** is on (for example) usefulness in wire layout. + ** + ***********************************************************************/ + +/*---------------------------------------------------------- + * String to path / path to string conversion. + * External format: + * "((xcoord, ycoord),... )" + * "[(xcoord, ycoord),... ]" + * "(xcoord, ycoord),... " + * "[xcoord, ycoord,... ]" + * Also support older format: + * "(closed, npts, xcoord, ycoord,... )" + *---------------------------------------------------------*/ + +Datum +path_area(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P(0); + float8 area = 0.0; + int i, + j; + + if (!path->closed) + PG_RETURN_NULL(); + + for (i = 0; i < path->npts; i++) + { + j = (i + 1) % path->npts; + area = float8_pl(area, float8_mul(path->p[i].x, path->p[j].y)); + area = float8_mi(area, float8_mul(path->p[i].y, path->p[j].x)); + } + + PG_RETURN_FLOAT8(float8_div(fabs(area), 2.0)); +} + + +Datum +path_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + PATH *path; + bool isopen; + char *s; + int npts; + int size; + int base_size; + int depth = 0; + + if ((npts = pair_count(str, ',')) <= 0) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + "path", str))); + + s = str; + while (isspace((unsigned char) *s)) + s++; + + /* skip single leading paren */ + if ((*s == LDELIM) && (strrchr(s, LDELIM) == s)) + { + s++; + depth++; + } + + base_size = sizeof(path->p[0]) * npts; + size = offsetof(PATH, p) + base_size; + + /* Check for integer overflow */ + if (base_size / npts != sizeof(path->p[0]) || size <= base_size) + ereport(ERROR, + (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED), + errmsg("too many points requested"))); + + path = (PATH *) palloc(size); + + SET_VARSIZE(path, size); + path->npts = npts; + + path_decode(s, true, npts, &(path->p[0]), &isopen, &s, "path", str); + + if (depth >= 1) + { + if (*s++ != RDELIM) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + "path", str))); + while (isspace((unsigned char) *s)) + s++; + } + if (*s != '\0') + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + "path", str))); + + path->closed = (!isopen); + /* prevent instability in unused pad bytes */ + path->dummy = 0; + + PG_RETURN_PATH_P(path); +} + + +Datum +path_out(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P(0); + + PG_RETURN_CSTRING(path_encode(path->closed ? PATH_CLOSED : PATH_OPEN, path->npts, path->p)); +} + +/* + * path_recv - converts external binary format to path + * + * External representation is closed flag (a boolean byte), int32 number + * of points, and the points. + */ +Datum +path_recv(PG_FUNCTION_ARGS) +{ + StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); + PATH *path; + int closed; + int32 npts; + int32 i; + int size; + + closed = pq_getmsgbyte(buf); + npts = pq_getmsgint(buf, sizeof(int32)); + if (npts <= 0 || npts >= (int32) ((INT_MAX - offsetof(PATH, p)) / sizeof(Point))) + ereport(ERROR, + (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION), + errmsg("invalid number of points in external \"path\" value"))); + + size = offsetof(PATH, p) + sizeof(path->p[0]) * npts; + path = (PATH *) palloc(size); + + SET_VARSIZE(path, size); + path->npts = npts; + path->closed = (closed ? 1 : 0); + /* prevent instability in unused pad bytes */ + path->dummy = 0; + + for (i = 0; i < npts; i++) + { + path->p[i].x = pq_getmsgfloat8(buf); + path->p[i].y = pq_getmsgfloat8(buf); + } + + PG_RETURN_PATH_P(path); +} + +/* + * path_send - converts path to binary format + */ +Datum +path_send(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P(0); + StringInfoData buf; + int32 i; + + pq_begintypsend(&buf); + pq_sendbyte(&buf, path->closed ? 1 : 0); + pq_sendint32(&buf, path->npts); + for (i = 0; i < path->npts; i++) + { + pq_sendfloat8(&buf, path->p[i].x); + pq_sendfloat8(&buf, path->p[i].y); + } + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); +} + + +/*---------------------------------------------------------- + * Relational operators. + * These are based on the path cardinality, + * as stupid as that sounds. + * + * Better relops and access methods coming soon. + *---------------------------------------------------------*/ + +Datum +path_n_lt(PG_FUNCTION_ARGS) +{ + PATH *p1 = PG_GETARG_PATH_P(0); + PATH *p2 = PG_GETARG_PATH_P(1); + + PG_RETURN_BOOL(p1->npts < p2->npts); +} + +Datum +path_n_gt(PG_FUNCTION_ARGS) +{ + PATH *p1 = PG_GETARG_PATH_P(0); + PATH *p2 = PG_GETARG_PATH_P(1); + + PG_RETURN_BOOL(p1->npts > p2->npts); +} + +Datum +path_n_eq(PG_FUNCTION_ARGS) +{ + PATH *p1 = PG_GETARG_PATH_P(0); + PATH *p2 = PG_GETARG_PATH_P(1); + + PG_RETURN_BOOL(p1->npts == p2->npts); +} + +Datum +path_n_le(PG_FUNCTION_ARGS) +{ + PATH *p1 = PG_GETARG_PATH_P(0); + PATH *p2 = PG_GETARG_PATH_P(1); + + PG_RETURN_BOOL(p1->npts <= p2->npts); +} + +Datum +path_n_ge(PG_FUNCTION_ARGS) +{ + PATH *p1 = PG_GETARG_PATH_P(0); + PATH *p2 = PG_GETARG_PATH_P(1); + + PG_RETURN_BOOL(p1->npts >= p2->npts); +} + +/*---------------------------------------------------------- + * Conversion operators. + *---------------------------------------------------------*/ + +Datum +path_isclosed(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P(0); + + PG_RETURN_BOOL(path->closed); +} + +Datum +path_isopen(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P(0); + + PG_RETURN_BOOL(!path->closed); +} + +Datum +path_npoints(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P(0); + + PG_RETURN_INT32(path->npts); +} + + +Datum +path_close(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P_COPY(0); + + path->closed = true; + + PG_RETURN_PATH_P(path); +} + +Datum +path_open(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P_COPY(0); + + path->closed = false; + + PG_RETURN_PATH_P(path); +} + + +/* path_inter - + * Does p1 intersect p2 at any point? + * Use bounding boxes for a quick (O(n)) check, then do a + * O(n^2) iterative edge check. + */ +Datum +path_inter(PG_FUNCTION_ARGS) +{ + PATH *p1 = PG_GETARG_PATH_P(0); + PATH *p2 = PG_GETARG_PATH_P(1); + BOX b1, + b2; + int i, + j; + LSEG seg1, + seg2; + + Assert(p1->npts > 0 && p2->npts > 0); + + b1.high.x = b1.low.x = p1->p[0].x; + b1.high.y = b1.low.y = p1->p[0].y; + for (i = 1; i < p1->npts; i++) + { + b1.high.x = float8_max(p1->p[i].x, b1.high.x); + b1.high.y = float8_max(p1->p[i].y, b1.high.y); + b1.low.x = float8_min(p1->p[i].x, b1.low.x); + b1.low.y = float8_min(p1->p[i].y, b1.low.y); + } + b2.high.x = b2.low.x = p2->p[0].x; + b2.high.y = b2.low.y = p2->p[0].y; + for (i = 1; i < p2->npts; i++) + { + b2.high.x = float8_max(p2->p[i].x, b2.high.x); + b2.high.y = float8_max(p2->p[i].y, b2.high.y); + b2.low.x = float8_min(p2->p[i].x, b2.low.x); + b2.low.y = float8_min(p2->p[i].y, b2.low.y); + } + if (!box_ov(&b1, &b2)) + PG_RETURN_BOOL(false); + + /* pairwise check lseg intersections */ + for (i = 0; i < p1->npts; i++) + { + int iprev; + + if (i > 0) + iprev = i - 1; + else + { + if (!p1->closed) + continue; + iprev = p1->npts - 1; /* include the closure segment */ + } + + for (j = 0; j < p2->npts; j++) + { + int jprev; + + if (j > 0) + jprev = j - 1; + else + { + if (!p2->closed) + continue; + jprev = p2->npts - 1; /* include the closure segment */ + } + + statlseg_construct(&seg1, &p1->p[iprev], &p1->p[i]); + statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]); + if (lseg_interpt_lseg(NULL, &seg1, &seg2)) + PG_RETURN_BOOL(true); + } + } + + /* if we dropped through, no two segs intersected */ + PG_RETURN_BOOL(false); +} + +/* path_distance() + * This essentially does a cartesian product of the lsegs in the + * two paths, and finds the min distance between any two lsegs + */ +Datum +path_distance(PG_FUNCTION_ARGS) +{ + PATH *p1 = PG_GETARG_PATH_P(0); + PATH *p2 = PG_GETARG_PATH_P(1); + float8 min = 0.0; /* initialize to keep compiler quiet */ + bool have_min = false; + float8 tmp; + int i, + j; + LSEG seg1, + seg2; + + for (i = 0; i < p1->npts; i++) + { + int iprev; + + if (i > 0) + iprev = i - 1; + else + { + if (!p1->closed) + continue; + iprev = p1->npts - 1; /* include the closure segment */ + } + + for (j = 0; j < p2->npts; j++) + { + int jprev; + + if (j > 0) + jprev = j - 1; + else + { + if (!p2->closed) + continue; + jprev = p2->npts - 1; /* include the closure segment */ + } + + statlseg_construct(&seg1, &p1->p[iprev], &p1->p[i]); + statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]); + + tmp = lseg_closept_lseg(NULL, &seg1, &seg2); + if (!have_min || float8_lt(tmp, min)) + { + min = tmp; + have_min = true; + } + } + } + + if (!have_min) + PG_RETURN_NULL(); + + PG_RETURN_FLOAT8(min); +} + + +/*---------------------------------------------------------- + * "Arithmetic" operations. + *---------------------------------------------------------*/ + +Datum +path_length(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P(0); + float8 result = 0.0; + int i; + + for (i = 0; i < path->npts; i++) + { + int iprev; + + if (i > 0) + iprev = i - 1; + else + { + if (!path->closed) + continue; + iprev = path->npts - 1; /* include the closure segment */ + } + + result = float8_pl(result, point_dt(&path->p[iprev], &path->p[i])); + } + + PG_RETURN_FLOAT8(result); +} + +/*********************************************************************** + ** + ** Routines for 2D points. + ** + ***********************************************************************/ + +/*---------------------------------------------------------- + * String to point, point to string conversion. + * External format: + * "(x,y)" + * "x,y" + *---------------------------------------------------------*/ + +Datum +point_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + Point *point = (Point *) palloc(sizeof(Point)); + + pair_decode(str, &point->x, &point->y, NULL, "point", str); + PG_RETURN_POINT_P(point); +} + +Datum +point_out(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + + PG_RETURN_CSTRING(path_encode(PATH_NONE, 1, pt)); +} + +/* + * point_recv - converts external binary format to point + */ +Datum +point_recv(PG_FUNCTION_ARGS) +{ + StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); + Point *point; + + point = (Point *) palloc(sizeof(Point)); + point->x = pq_getmsgfloat8(buf); + point->y = pq_getmsgfloat8(buf); + PG_RETURN_POINT_P(point); +} + +/* + * point_send - converts point to binary format + */ +Datum +point_send(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + StringInfoData buf; + + pq_begintypsend(&buf); + pq_sendfloat8(&buf, pt->x); + pq_sendfloat8(&buf, pt->y); + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); +} + + +/* + * Initialize a point + */ +static inline void +point_construct(Point *result, float8 x, float8 y) +{ + result->x = x; + result->y = y; +} + + +/*---------------------------------------------------------- + * Relational operators for Points. + * Since we do have a sense of coordinates being + * "equal" to a given accuracy (point_vert, point_horiz), + * the other ops must preserve that sense. This means + * that results may, strictly speaking, be a lie (unless + * EPSILON = 0.0). + *---------------------------------------------------------*/ + +Datum +point_left(PG_FUNCTION_ARGS) +{ + Point *pt1 = PG_GETARG_POINT_P(0); + Point *pt2 = PG_GETARG_POINT_P(1); + + PG_RETURN_BOOL(FPlt(pt1->x, pt2->x)); +} + +Datum +point_right(PG_FUNCTION_ARGS) +{ + Point *pt1 = PG_GETARG_POINT_P(0); + Point *pt2 = PG_GETARG_POINT_P(1); + + PG_RETURN_BOOL(FPgt(pt1->x, pt2->x)); +} + +Datum +point_above(PG_FUNCTION_ARGS) +{ + Point *pt1 = PG_GETARG_POINT_P(0); + Point *pt2 = PG_GETARG_POINT_P(1); + + PG_RETURN_BOOL(FPgt(pt1->y, pt2->y)); +} + +Datum +point_below(PG_FUNCTION_ARGS) +{ + Point *pt1 = PG_GETARG_POINT_P(0); + Point *pt2 = PG_GETARG_POINT_P(1); + + PG_RETURN_BOOL(FPlt(pt1->y, pt2->y)); +} + +Datum +point_vert(PG_FUNCTION_ARGS) +{ + Point *pt1 = PG_GETARG_POINT_P(0); + Point *pt2 = PG_GETARG_POINT_P(1); + + PG_RETURN_BOOL(FPeq(pt1->x, pt2->x)); +} + +Datum +point_horiz(PG_FUNCTION_ARGS) +{ + Point *pt1 = PG_GETARG_POINT_P(0); + Point *pt2 = PG_GETARG_POINT_P(1); + + PG_RETURN_BOOL(FPeq(pt1->y, pt2->y)); +} + +Datum +point_eq(PG_FUNCTION_ARGS) +{ + Point *pt1 = PG_GETARG_POINT_P(0); + Point *pt2 = PG_GETARG_POINT_P(1); + + PG_RETURN_BOOL(point_eq_point(pt1, pt2)); +} + +Datum +point_ne(PG_FUNCTION_ARGS) +{ + Point *pt1 = PG_GETARG_POINT_P(0); + Point *pt2 = PG_GETARG_POINT_P(1); + + PG_RETURN_BOOL(!point_eq_point(pt1, pt2)); +} + + +/* + * Check whether the two points are the same + */ +static inline bool +point_eq_point(Point *pt1, Point *pt2) +{ + /* If any NaNs are involved, insist on exact equality */ + if (unlikely(isnan(pt1->x) || isnan(pt1->y) || + isnan(pt2->x) || isnan(pt2->y))) + return (float8_eq(pt1->x, pt2->x) && float8_eq(pt1->y, pt2->y)); + + return (FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y)); +} + + +/*---------------------------------------------------------- + * "Arithmetic" operators on points. + *---------------------------------------------------------*/ + +Datum +point_distance(PG_FUNCTION_ARGS) +{ + Point *pt1 = PG_GETARG_POINT_P(0); + Point *pt2 = PG_GETARG_POINT_P(1); + + PG_RETURN_FLOAT8(point_dt(pt1, pt2)); +} + +static inline float8 +point_dt(Point *pt1, Point *pt2) +{ + return HYPOT(float8_mi(pt1->x, pt2->x), float8_mi(pt1->y, pt2->y)); +} + +Datum +point_slope(PG_FUNCTION_ARGS) +{ + Point *pt1 = PG_GETARG_POINT_P(0); + Point *pt2 = PG_GETARG_POINT_P(1); + + PG_RETURN_FLOAT8(point_sl(pt1, pt2)); +} + + +/* + * Return slope of two points + * + * Note that this function returns Inf when the points are the same. + */ +static inline float8 +point_sl(Point *pt1, Point *pt2) +{ + if (FPeq(pt1->x, pt2->x)) + return get_float8_infinity(); + if (FPeq(pt1->y, pt2->y)) + return 0.0; + return float8_div(float8_mi(pt1->y, pt2->y), float8_mi(pt1->x, pt2->x)); +} + + +/* + * Return inverse slope of two points + * + * Note that this function returns 0.0 when the points are the same. + */ +static inline float8 +point_invsl(Point *pt1, Point *pt2) +{ + if (FPeq(pt1->x, pt2->x)) + return 0.0; + if (FPeq(pt1->y, pt2->y)) + return get_float8_infinity(); + return float8_div(float8_mi(pt1->x, pt2->x), float8_mi(pt2->y, pt1->y)); +} + + +/*********************************************************************** + ** + ** Routines for 2D line segments. + ** + ***********************************************************************/ + +/*---------------------------------------------------------- + * String to lseg, lseg to string conversion. + * External forms: "[(x1, y1), (x2, y2)]" + * "(x1, y1), (x2, y2)" + * "x1, y1, x2, y2" + * closed form ok "((x1, y1), (x2, y2))" + * (old form) "(x1, y1, x2, y2)" + *---------------------------------------------------------*/ + +Datum +lseg_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + LSEG *lseg = (LSEG *) palloc(sizeof(LSEG)); + bool isopen; + + path_decode(str, true, 2, &lseg->p[0], &isopen, NULL, "lseg", str); + PG_RETURN_LSEG_P(lseg); +} + + +Datum +lseg_out(PG_FUNCTION_ARGS) +{ + LSEG *ls = PG_GETARG_LSEG_P(0); + + PG_RETURN_CSTRING(path_encode(PATH_OPEN, 2, &ls->p[0])); +} + +/* + * lseg_recv - converts external binary format to lseg + */ +Datum +lseg_recv(PG_FUNCTION_ARGS) +{ + StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); + LSEG *lseg; + + lseg = (LSEG *) palloc(sizeof(LSEG)); + + lseg->p[0].x = pq_getmsgfloat8(buf); + lseg->p[0].y = pq_getmsgfloat8(buf); + lseg->p[1].x = pq_getmsgfloat8(buf); + lseg->p[1].y = pq_getmsgfloat8(buf); + + PG_RETURN_LSEG_P(lseg); +} + +/* + * lseg_send - converts lseg to binary format + */ +Datum +lseg_send(PG_FUNCTION_ARGS) +{ + LSEG *ls = PG_GETARG_LSEG_P(0); + StringInfoData buf; + + pq_begintypsend(&buf); + pq_sendfloat8(&buf, ls->p[0].x); + pq_sendfloat8(&buf, ls->p[0].y); + pq_sendfloat8(&buf, ls->p[1].x); + pq_sendfloat8(&buf, ls->p[1].y); + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); +} + + +/* lseg_construct - + * form a LSEG from two Points. + */ +Datum +lseg_construct(PG_FUNCTION_ARGS) +{ + Point *pt1 = PG_GETARG_POINT_P(0); + Point *pt2 = PG_GETARG_POINT_P(1); + LSEG *result = (LSEG *) palloc(sizeof(LSEG)); + + statlseg_construct(result, pt1, pt2); + + PG_RETURN_LSEG_P(result); +} + +/* like lseg_construct, but assume space already allocated */ +static inline void +statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2) +{ + lseg->p[0].x = pt1->x; + lseg->p[0].y = pt1->y; + lseg->p[1].x = pt2->x; + lseg->p[1].y = pt2->y; +} + + +/* + * Return slope of the line segment + */ +static inline float8 +lseg_sl(LSEG *lseg) +{ + return point_sl(&lseg->p[0], &lseg->p[1]); +} + + +/* + * Return inverse slope of the line segment + */ +static inline float8 +lseg_invsl(LSEG *lseg) +{ + return point_invsl(&lseg->p[0], &lseg->p[1]); +} + + +Datum +lseg_length(PG_FUNCTION_ARGS) +{ + LSEG *lseg = PG_GETARG_LSEG_P(0); + + PG_RETURN_FLOAT8(point_dt(&lseg->p[0], &lseg->p[1])); +} + +/*---------------------------------------------------------- + * Relative position routines. + *---------------------------------------------------------*/ + +/* + ** find intersection of the two lines, and see if it falls on + ** both segments. + */ +Datum +lseg_intersect(PG_FUNCTION_ARGS) +{ + LSEG *l1 = PG_GETARG_LSEG_P(0); + LSEG *l2 = PG_GETARG_LSEG_P(1); + + PG_RETURN_BOOL(lseg_interpt_lseg(NULL, l1, l2)); +} + + +Datum +lseg_parallel(PG_FUNCTION_ARGS) +{ + LSEG *l1 = PG_GETARG_LSEG_P(0); + LSEG *l2 = PG_GETARG_LSEG_P(1); + + PG_RETURN_BOOL(FPeq(lseg_sl(l1), lseg_sl(l2))); +} + +/* + * Determine if two line segments are perpendicular. + */ +Datum +lseg_perp(PG_FUNCTION_ARGS) +{ + LSEG *l1 = PG_GETARG_LSEG_P(0); + LSEG *l2 = PG_GETARG_LSEG_P(1); + + PG_RETURN_BOOL(FPeq(lseg_sl(l1), lseg_invsl(l2))); +} + +Datum +lseg_vertical(PG_FUNCTION_ARGS) +{ + LSEG *lseg = PG_GETARG_LSEG_P(0); + + PG_RETURN_BOOL(FPeq(lseg->p[0].x, lseg->p[1].x)); +} + +Datum +lseg_horizontal(PG_FUNCTION_ARGS) +{ + LSEG *lseg = PG_GETARG_LSEG_P(0); + + PG_RETURN_BOOL(FPeq(lseg->p[0].y, lseg->p[1].y)); +} + + +Datum +lseg_eq(PG_FUNCTION_ARGS) +{ + LSEG *l1 = PG_GETARG_LSEG_P(0); + LSEG *l2 = PG_GETARG_LSEG_P(1); + + PG_RETURN_BOOL(point_eq_point(&l1->p[0], &l2->p[0]) && + point_eq_point(&l1->p[1], &l2->p[1])); +} + +Datum +lseg_ne(PG_FUNCTION_ARGS) +{ + LSEG *l1 = PG_GETARG_LSEG_P(0); + LSEG *l2 = PG_GETARG_LSEG_P(1); + + PG_RETURN_BOOL(!point_eq_point(&l1->p[0], &l2->p[0]) || + !point_eq_point(&l1->p[1], &l2->p[1])); +} + +Datum +lseg_lt(PG_FUNCTION_ARGS) +{ + LSEG *l1 = PG_GETARG_LSEG_P(0); + LSEG *l2 = PG_GETARG_LSEG_P(1); + + PG_RETURN_BOOL(FPlt(point_dt(&l1->p[0], &l1->p[1]), + point_dt(&l2->p[0], &l2->p[1]))); +} + +Datum +lseg_le(PG_FUNCTION_ARGS) +{ + LSEG *l1 = PG_GETARG_LSEG_P(0); + LSEG *l2 = PG_GETARG_LSEG_P(1); + + PG_RETURN_BOOL(FPle(point_dt(&l1->p[0], &l1->p[1]), + point_dt(&l2->p[0], &l2->p[1]))); +} + +Datum +lseg_gt(PG_FUNCTION_ARGS) +{ + LSEG *l1 = PG_GETARG_LSEG_P(0); + LSEG *l2 = PG_GETARG_LSEG_P(1); + + PG_RETURN_BOOL(FPgt(point_dt(&l1->p[0], &l1->p[1]), + point_dt(&l2->p[0], &l2->p[1]))); +} + +Datum +lseg_ge(PG_FUNCTION_ARGS) +{ + LSEG *l1 = PG_GETARG_LSEG_P(0); + LSEG *l2 = PG_GETARG_LSEG_P(1); + + PG_RETURN_BOOL(FPge(point_dt(&l1->p[0], &l1->p[1]), + point_dt(&l2->p[0], &l2->p[1]))); +} + + +/*---------------------------------------------------------- + * Line arithmetic routines. + *---------------------------------------------------------*/ + +/* lseg_distance - + * If two segments don't intersect, then the closest + * point will be from one of the endpoints to the other + * segment. + */ +Datum +lseg_distance(PG_FUNCTION_ARGS) +{ + LSEG *l1 = PG_GETARG_LSEG_P(0); + LSEG *l2 = PG_GETARG_LSEG_P(1); + + PG_RETURN_FLOAT8(lseg_closept_lseg(NULL, l1, l2)); +} + + +Datum +lseg_center(PG_FUNCTION_ARGS) +{ + LSEG *lseg = PG_GETARG_LSEG_P(0); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + result->x = float8_div(float8_pl(lseg->p[0].x, lseg->p[1].x), 2.0); + result->y = float8_div(float8_pl(lseg->p[0].y, lseg->p[1].y), 2.0); + + PG_RETURN_POINT_P(result); +} + + +/* + * Return whether the two segments intersect. If *result is not NULL, + * it is set to the intersection point. + * + * This function is almost perfectly symmetric, even though it doesn't look + * like it. See lseg_interpt_line() for the other half of it. + */ +static bool +lseg_interpt_lseg(Point *result, LSEG *l1, LSEG *l2) +{ + Point interpt; + LINE tmp; + + line_construct(&tmp, &l2->p[0], lseg_sl(l2)); + if (!lseg_interpt_line(&interpt, l1, &tmp)) + return false; + + /* + * If the line intersection point isn't within l2, there is no valid + * segment intersection point at all. + */ + if (!lseg_contain_point(l2, &interpt)) + return false; + + if (result != NULL) + *result = interpt; + + return true; +} + +Datum +lseg_interpt(PG_FUNCTION_ARGS) +{ + LSEG *l1 = PG_GETARG_LSEG_P(0); + LSEG *l2 = PG_GETARG_LSEG_P(1); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + if (!lseg_interpt_lseg(result, l1, l2)) + PG_RETURN_NULL(); + PG_RETURN_POINT_P(result); +} + +/*********************************************************************** + ** + ** Routines for position comparisons of differently-typed + ** 2D objects. + ** + ***********************************************************************/ + +/*--------------------------------------------------------------------- + * dist_ + * Minimum distance from one object to another. + *-------------------------------------------------------------------*/ + +/* + * Distance from a point to a line + */ +Datum +dist_pl(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + LINE *line = PG_GETARG_LINE_P(1); + + PG_RETURN_FLOAT8(line_closept_point(NULL, line, pt)); +} + +/* + * Distance from a line to a point + */ +Datum +dist_lp(PG_FUNCTION_ARGS) +{ + LINE *line = PG_GETARG_LINE_P(0); + Point *pt = PG_GETARG_POINT_P(1); + + PG_RETURN_FLOAT8(line_closept_point(NULL, line, pt)); +} + +/* + * Distance from a point to a lseg + */ +Datum +dist_ps(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + LSEG *lseg = PG_GETARG_LSEG_P(1); + + PG_RETURN_FLOAT8(lseg_closept_point(NULL, lseg, pt)); +} + +/* + * Distance from a lseg to a point + */ +Datum +dist_sp(PG_FUNCTION_ARGS) +{ + LSEG *lseg = PG_GETARG_LSEG_P(0); + Point *pt = PG_GETARG_POINT_P(1); + + PG_RETURN_FLOAT8(lseg_closept_point(NULL, lseg, pt)); +} + +static float8 +dist_ppath_internal(Point *pt, PATH *path) +{ + float8 result = 0.0; /* keep compiler quiet */ + bool have_min = false; + float8 tmp; + int i; + LSEG lseg; + + Assert(path->npts > 0); + + /* + * The distance from a point to a path is the smallest distance from the + * point to any of its constituent segments. + */ + for (i = 0; i < path->npts; i++) + { + int iprev; + + if (i > 0) + iprev = i - 1; + else + { + if (!path->closed) + continue; + iprev = path->npts - 1; /* Include the closure segment */ + } + + statlseg_construct(&lseg, &path->p[iprev], &path->p[i]); + tmp = lseg_closept_point(NULL, &lseg, pt); + if (!have_min || float8_lt(tmp, result)) + { + result = tmp; + have_min = true; + } + } + + return result; +} + +/* + * Distance from a point to a path + */ +Datum +dist_ppath(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + PATH *path = PG_GETARG_PATH_P(1); + + PG_RETURN_FLOAT8(dist_ppath_internal(pt, path)); +} + +/* + * Distance from a path to a point + */ +Datum +dist_pathp(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P(0); + Point *pt = PG_GETARG_POINT_P(1); + + PG_RETURN_FLOAT8(dist_ppath_internal(pt, path)); +} + +/* + * Distance from a point to a box + */ +Datum +dist_pb(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + BOX *box = PG_GETARG_BOX_P(1); + + PG_RETURN_FLOAT8(box_closept_point(NULL, box, pt)); +} + +/* + * Distance from a box to a point + */ +Datum +dist_bp(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + Point *pt = PG_GETARG_POINT_P(1); + + PG_RETURN_FLOAT8(box_closept_point(NULL, box, pt)); +} + +/* + * Distance from a lseg to a line + */ +Datum +dist_sl(PG_FUNCTION_ARGS) +{ + LSEG *lseg = PG_GETARG_LSEG_P(0); + LINE *line = PG_GETARG_LINE_P(1); + + PG_RETURN_FLOAT8(lseg_closept_line(NULL, lseg, line)); +} + +/* + * Distance from a line to a lseg + */ +Datum +dist_ls(PG_FUNCTION_ARGS) +{ + LINE *line = PG_GETARG_LINE_P(0); + LSEG *lseg = PG_GETARG_LSEG_P(1); + + PG_RETURN_FLOAT8(lseg_closept_line(NULL, lseg, line)); +} + +/* + * Distance from a lseg to a box + */ +Datum +dist_sb(PG_FUNCTION_ARGS) +{ + LSEG *lseg = PG_GETARG_LSEG_P(0); + BOX *box = PG_GETARG_BOX_P(1); + + PG_RETURN_FLOAT8(box_closept_lseg(NULL, box, lseg)); +} + +/* + * Distance from a box to a lseg + */ +Datum +dist_bs(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + LSEG *lseg = PG_GETARG_LSEG_P(1); + + PG_RETURN_FLOAT8(box_closept_lseg(NULL, box, lseg)); +} + +/* + * Distance from a line to a box + */ +Datum +dist_lb(PG_FUNCTION_ARGS) +{ +#ifdef NOT_USED + LINE *line = PG_GETARG_LINE_P(0); + BOX *box = PG_GETARG_BOX_P(1); +#endif + + /* need to think about this one for a while */ + ereport(ERROR, + (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), + errmsg("function \"dist_lb\" not implemented"))); + + PG_RETURN_NULL(); +} + +/* + * Distance from a box to a line + */ +Datum +dist_bl(PG_FUNCTION_ARGS) +{ +#ifdef NOT_USED + BOX *box = PG_GETARG_BOX_P(0); + LINE *line = PG_GETARG_LINE_P(1); +#endif + + /* need to think about this one for a while */ + ereport(ERROR, + (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), + errmsg("function \"dist_bl\" not implemented"))); + + PG_RETURN_NULL(); +} + +static float8 +dist_cpoly_internal(CIRCLE *circle, POLYGON *poly) +{ + float8 result; + + /* calculate distance to center, and subtract radius */ + result = float8_mi(dist_ppoly_internal(&circle->center, poly), + circle->radius); + if (result < 0.0) + result = 0.0; + + return result; +} + +/* + * Distance from a circle to a polygon + */ +Datum +dist_cpoly(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + POLYGON *poly = PG_GETARG_POLYGON_P(1); + + PG_RETURN_FLOAT8(dist_cpoly_internal(circle, poly)); +} + +/* + * Distance from a polygon to a circle + */ +Datum +dist_polyc(PG_FUNCTION_ARGS) +{ + POLYGON *poly = PG_GETARG_POLYGON_P(0); + CIRCLE *circle = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_FLOAT8(dist_cpoly_internal(circle, poly)); +} + +/* + * Distance from a point to a polygon + */ +Datum +dist_ppoly(PG_FUNCTION_ARGS) +{ + Point *point = PG_GETARG_POINT_P(0); + POLYGON *poly = PG_GETARG_POLYGON_P(1); + + PG_RETURN_FLOAT8(dist_ppoly_internal(point, poly)); +} + +Datum +dist_polyp(PG_FUNCTION_ARGS) +{ + POLYGON *poly = PG_GETARG_POLYGON_P(0); + Point *point = PG_GETARG_POINT_P(1); + + PG_RETURN_FLOAT8(dist_ppoly_internal(point, poly)); +} + +static float8 +dist_ppoly_internal(Point *pt, POLYGON *poly) +{ + float8 result; + float8 d; + int i; + LSEG seg; + + if (point_inside(pt, poly->npts, poly->p) != 0) + return 0.0; + + /* initialize distance with segment between first and last points */ + seg.p[0].x = poly->p[0].x; + seg.p[0].y = poly->p[0].y; + seg.p[1].x = poly->p[poly->npts - 1].x; + seg.p[1].y = poly->p[poly->npts - 1].y; + result = lseg_closept_point(NULL, &seg, pt); + + /* check distances for other segments */ + for (i = 0; i < poly->npts - 1; i++) + { + seg.p[0].x = poly->p[i].x; + seg.p[0].y = poly->p[i].y; + seg.p[1].x = poly->p[i + 1].x; + seg.p[1].y = poly->p[i + 1].y; + d = lseg_closept_point(NULL, &seg, pt); + if (float8_lt(d, result)) + result = d; + } + + return result; +} + + +/*--------------------------------------------------------------------- + * interpt_ + * Intersection point of objects. + * We choose to ignore the "point" of intersection between + * lines and boxes, since there are typically two. + *-------------------------------------------------------------------*/ + +/* + * Return whether the line segment intersect with the line. If *result is not + * NULL, it is set to the intersection point. + */ +static bool +lseg_interpt_line(Point *result, LSEG *lseg, LINE *line) +{ + Point interpt; + LINE tmp; + + /* + * First, we promote the line segment to a line, because we know how to + * find the intersection point of two lines. If they don't have an + * intersection point, we are done. + */ + line_construct(&tmp, &lseg->p[0], lseg_sl(lseg)); + if (!line_interpt_line(&interpt, &tmp, line)) + return false; + + /* + * Then, we check whether the intersection point is actually on the line + * segment. + */ + if (!lseg_contain_point(lseg, &interpt)) + return false; + if (result != NULL) + { + /* + * If there is an intersection, then check explicitly for matching + * endpoints since there may be rounding effects with annoying LSB + * residue. + */ + if (point_eq_point(&lseg->p[0], &interpt)) + *result = lseg->p[0]; + else if (point_eq_point(&lseg->p[1], &interpt)) + *result = lseg->p[1]; + else + *result = interpt; + } + + return true; +} + +/*--------------------------------------------------------------------- + * close_ + * Point of closest proximity between objects. + *-------------------------------------------------------------------*/ + +/* + * If *result is not NULL, it is set to the intersection point of a + * perpendicular of the line through the point. Returns the distance + * of those two points. + */ +static float8 +line_closept_point(Point *result, LINE *line, Point *point) +{ + Point closept; + LINE tmp; + + /* + * We drop a perpendicular to find the intersection point. Ordinarily we + * should always find it, but that can fail in the presence of NaN + * coordinates, and perhaps even from simple roundoff issues. + */ + line_construct(&tmp, point, line_invsl(line)); + if (!line_interpt_line(&closept, &tmp, line)) + { + if (result != NULL) + *result = *point; + + return get_float8_nan(); + } + + if (result != NULL) + *result = closept; + + return point_dt(&closept, point); +} + +Datum +close_pl(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + LINE *line = PG_GETARG_LINE_P(1); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + if (isnan(line_closept_point(result, line, pt))) + PG_RETURN_NULL(); + + PG_RETURN_POINT_P(result); +} + + +/* + * Closest point on line segment to specified point. + * + * If *result is not NULL, set it to the closest point on the line segment + * to the point. Returns the distance of the two points. + */ +static float8 +lseg_closept_point(Point *result, LSEG *lseg, Point *pt) +{ + Point closept; + LINE tmp; + + /* + * To find the closest point, we draw a perpendicular line from the point + * to the line segment. + */ + line_construct(&tmp, pt, point_invsl(&lseg->p[0], &lseg->p[1])); + lseg_closept_line(&closept, lseg, &tmp); + + if (result != NULL) + *result = closept; + + return point_dt(&closept, pt); +} + +Datum +close_ps(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + LSEG *lseg = PG_GETARG_LSEG_P(1); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + if (isnan(lseg_closept_point(result, lseg, pt))) + PG_RETURN_NULL(); + + PG_RETURN_POINT_P(result); +} + + +/* + * Closest point on line segment to line segment + */ +static float8 +lseg_closept_lseg(Point *result, LSEG *on_lseg, LSEG *to_lseg) +{ + Point point; + float8 dist, + d; + + /* First, we handle the case when the line segments are intersecting. */ + if (lseg_interpt_lseg(result, on_lseg, to_lseg)) + return 0.0; + + /* + * Then, we find the closest points from the endpoints of the second line + * segment, and keep the closest one. + */ + dist = lseg_closept_point(result, on_lseg, &to_lseg->p[0]); + d = lseg_closept_point(&point, on_lseg, &to_lseg->p[1]); + if (float8_lt(d, dist)) + { + dist = d; + if (result != NULL) + *result = point; + } + + /* The closest point can still be one of the endpoints, so we test them. */ + d = lseg_closept_point(NULL, to_lseg, &on_lseg->p[0]); + if (float8_lt(d, dist)) + { + dist = d; + if (result != NULL) + *result = on_lseg->p[0]; + } + d = lseg_closept_point(NULL, to_lseg, &on_lseg->p[1]); + if (float8_lt(d, dist)) + { + dist = d; + if (result != NULL) + *result = on_lseg->p[1]; + } + + return dist; +} + +Datum +close_lseg(PG_FUNCTION_ARGS) +{ + LSEG *l1 = PG_GETARG_LSEG_P(0); + LSEG *l2 = PG_GETARG_LSEG_P(1); + Point *result; + + if (lseg_sl(l1) == lseg_sl(l2)) + PG_RETURN_NULL(); + + result = (Point *) palloc(sizeof(Point)); + + if (isnan(lseg_closept_lseg(result, l2, l1))) + PG_RETURN_NULL(); + + PG_RETURN_POINT_P(result); +} + + +/* + * Closest point on or in box to specified point. + * + * If *result is not NULL, set it to the closest point on the box to the + * given point, and return the distance of the two points. + */ +static float8 +box_closept_point(Point *result, BOX *box, Point *pt) +{ + float8 dist, + d; + Point point, + closept; + LSEG lseg; + + if (box_contain_point(box, pt)) + { + if (result != NULL) + *result = *pt; + + return 0.0; + } + + /* pairwise check lseg distances */ + point.x = box->low.x; + point.y = box->high.y; + statlseg_construct(&lseg, &box->low, &point); + dist = lseg_closept_point(result, &lseg, pt); + + statlseg_construct(&lseg, &box->high, &point); + d = lseg_closept_point(&closept, &lseg, pt); + if (float8_lt(d, dist)) + { + dist = d; + if (result != NULL) + *result = closept; + } + + point.x = box->high.x; + point.y = box->low.y; + statlseg_construct(&lseg, &box->low, &point); + d = lseg_closept_point(&closept, &lseg, pt); + if (float8_lt(d, dist)) + { + dist = d; + if (result != NULL) + *result = closept; + } + + statlseg_construct(&lseg, &box->high, &point); + d = lseg_closept_point(&closept, &lseg, pt); + if (float8_lt(d, dist)) + { + dist = d; + if (result != NULL) + *result = closept; + } + + return dist; +} + +Datum +close_pb(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + BOX *box = PG_GETARG_BOX_P(1); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + if (isnan(box_closept_point(result, box, pt))) + PG_RETURN_NULL(); + + PG_RETURN_POINT_P(result); +} + + +/* close_sl() + * Closest point on line to line segment. + * + * XXX THIS CODE IS WRONG + * The code is actually calculating the point on the line segment + * which is backwards from the routine naming convention. + * Copied code to new routine close_ls() but haven't fixed this one yet. + * - thomas 1998-01-31 + */ +Datum +close_sl(PG_FUNCTION_ARGS) +{ +#ifdef NOT_USED + LSEG *lseg = PG_GETARG_LSEG_P(0); + LINE *line = PG_GETARG_LINE_P(1); + Point *result; + float8 d1, + d2; + + result = (Point *) palloc(sizeof(Point)); + + if (lseg_interpt_line(result, lseg, line)) + PG_RETURN_POINT_P(result); + + d1 = line_closept_point(NULL, line, &lseg->p[0]); + d2 = line_closept_point(NULL, line, &lseg->p[1]); + if (float8_lt(d1, d2)) + *result = lseg->p[0]; + else + *result = lseg->p[1]; + + PG_RETURN_POINT_P(result); +#endif + + ereport(ERROR, + (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), + errmsg("function \"close_sl\" not implemented"))); + + PG_RETURN_NULL(); +} + +/* + * Closest point on line segment to line. + * + * Return the distance between the line and the closest point of the line + * segment to the line. If *result is not NULL, set it to that point. + * + * NOTE: When the lines are parallel, endpoints of one of the line segment + * are FPeq(), in presence of NaN or Infinite coordinates, or perhaps = + * even because of simple roundoff issues, there may not be a single closest + * point. We are likely to set the result to the second endpoint in these + * cases. + */ +static float8 +lseg_closept_line(Point *result, LSEG *lseg, LINE *line) +{ + float8 dist1, + dist2; + + if (lseg_interpt_line(result, lseg, line)) + return 0.0; + + dist1 = line_closept_point(NULL, line, &lseg->p[0]); + dist2 = line_closept_point(NULL, line, &lseg->p[1]); + + if (dist1 < dist2) + { + if (result != NULL) + *result = lseg->p[0]; + + return dist1; + } + else + { + if (result != NULL) + *result = lseg->p[1]; + + return dist2; + } +} + +Datum +close_ls(PG_FUNCTION_ARGS) +{ + LINE *line = PG_GETARG_LINE_P(0); + LSEG *lseg = PG_GETARG_LSEG_P(1); + Point *result; + + if (lseg_sl(lseg) == line_sl(line)) + PG_RETURN_NULL(); + + result = (Point *) palloc(sizeof(Point)); + + if (isnan(lseg_closept_line(result, lseg, line))) + PG_RETURN_NULL(); + + PG_RETURN_POINT_P(result); +} + + +/* + * Closest point on or in box to line segment. + * + * Returns the distance between the closest point on or in the box to + * the line segment. If *result is not NULL, it is set to that point. + */ +static float8 +box_closept_lseg(Point *result, BOX *box, LSEG *lseg) +{ + float8 dist, + d; + Point point, + closept; + LSEG bseg; + + if (box_interpt_lseg(result, box, lseg)) + return 0.0; + + /* pairwise check lseg distances */ + point.x = box->low.x; + point.y = box->high.y; + statlseg_construct(&bseg, &box->low, &point); + dist = lseg_closept_lseg(result, &bseg, lseg); + + statlseg_construct(&bseg, &box->high, &point); + d = lseg_closept_lseg(&closept, &bseg, lseg); + if (float8_lt(d, dist)) + { + dist = d; + if (result != NULL) + *result = closept; + } + + point.x = box->high.x; + point.y = box->low.y; + statlseg_construct(&bseg, &box->low, &point); + d = lseg_closept_lseg(&closept, &bseg, lseg); + if (float8_lt(d, dist)) + { + dist = d; + if (result != NULL) + *result = closept; + } + + statlseg_construct(&bseg, &box->high, &point); + d = lseg_closept_lseg(&closept, &bseg, lseg); + if (float8_lt(d, dist)) + { + dist = d; + if (result != NULL) + *result = closept; + } + + return dist; +} + +Datum +close_sb(PG_FUNCTION_ARGS) +{ + LSEG *lseg = PG_GETARG_LSEG_P(0); + BOX *box = PG_GETARG_BOX_P(1); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + if (isnan(box_closept_lseg(result, box, lseg))) + PG_RETURN_NULL(); + + PG_RETURN_POINT_P(result); +} + + +Datum +close_lb(PG_FUNCTION_ARGS) +{ +#ifdef NOT_USED + LINE *line = PG_GETARG_LINE_P(0); + BOX *box = PG_GETARG_BOX_P(1); +#endif + + /* think about this one for a while */ + ereport(ERROR, + (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), + errmsg("function \"close_lb\" not implemented"))); + + PG_RETURN_NULL(); +} + +/*--------------------------------------------------------------------- + * on_ + * Whether one object lies completely within another. + *-------------------------------------------------------------------*/ + +/* + * Does the point satisfy the equation? + */ +static bool +line_contain_point(LINE *line, Point *point) +{ + return FPzero(float8_pl(float8_pl(float8_mul(line->A, point->x), + float8_mul(line->B, point->y)), + line->C)); +} + +Datum +on_pl(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + LINE *line = PG_GETARG_LINE_P(1); + + PG_RETURN_BOOL(line_contain_point(line, pt)); +} + + +/* + * Determine colinearity by detecting a triangle inequality. + * This algorithm seems to behave nicely even with lsb residues - tgl 1997-07-09 + */ +static bool +lseg_contain_point(LSEG *lseg, Point *pt) +{ + return FPeq(point_dt(pt, &lseg->p[0]) + + point_dt(pt, &lseg->p[1]), + point_dt(&lseg->p[0], &lseg->p[1])); +} + +Datum +on_ps(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + LSEG *lseg = PG_GETARG_LSEG_P(1); + + PG_RETURN_BOOL(lseg_contain_point(lseg, pt)); +} + + +/* + * Check whether the point is in the box or on its border + */ +static bool +box_contain_point(BOX *box, Point *point) +{ + return box->high.x >= point->x && box->low.x <= point->x && + box->high.y >= point->y && box->low.y <= point->y; +} + +Datum +on_pb(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + BOX *box = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(box_contain_point(box, pt)); +} + +Datum +box_contain_pt(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + Point *pt = PG_GETARG_POINT_P(1); + + PG_RETURN_BOOL(box_contain_point(box, pt)); +} + +/* on_ppath - + * Whether a point lies within (on) a polyline. + * If open, we have to (groan) check each segment. + * (uses same algorithm as for point intersecting segment - tgl 1997-07-09) + * If closed, we use the old O(n) ray method for point-in-polygon. + * The ray is horizontal, from pt out to the right. + * Each segment that crosses the ray counts as an + * intersection; note that an endpoint or edge may touch + * but not cross. + * (we can do p-in-p in lg(n), but it takes preprocessing) + */ +Datum +on_ppath(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + PATH *path = PG_GETARG_PATH_P(1); + int i, + n; + float8 a, + b; + + /*-- OPEN --*/ + if (!path->closed) + { + n = path->npts - 1; + a = point_dt(pt, &path->p[0]); + for (i = 0; i < n; i++) + { + b = point_dt(pt, &path->p[i + 1]); + if (FPeq(float8_pl(a, b), point_dt(&path->p[i], &path->p[i + 1]))) + PG_RETURN_BOOL(true); + a = b; + } + PG_RETURN_BOOL(false); + } + + /*-- CLOSED --*/ + PG_RETURN_BOOL(point_inside(pt, path->npts, path->p) != 0); +} + + +/* + * Check whether the line segment is on the line or close enough + * + * It is, if both of its points are on the line or close enough. + */ +Datum +on_sl(PG_FUNCTION_ARGS) +{ + LSEG *lseg = PG_GETARG_LSEG_P(0); + LINE *line = PG_GETARG_LINE_P(1); + + PG_RETURN_BOOL(line_contain_point(line, &lseg->p[0]) && + line_contain_point(line, &lseg->p[1])); +} + + +/* + * Check whether the line segment is in the box or on its border + * + * It is, if both of its points are in the box or on its border. + */ +static bool +box_contain_lseg(BOX *box, LSEG *lseg) +{ + return box_contain_point(box, &lseg->p[0]) && + box_contain_point(box, &lseg->p[1]); +} + +Datum +on_sb(PG_FUNCTION_ARGS) +{ + LSEG *lseg = PG_GETARG_LSEG_P(0); + BOX *box = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(box_contain_lseg(box, lseg)); +} + +/*--------------------------------------------------------------------- + * inter_ + * Whether one object intersects another. + *-------------------------------------------------------------------*/ + +Datum +inter_sl(PG_FUNCTION_ARGS) +{ + LSEG *lseg = PG_GETARG_LSEG_P(0); + LINE *line = PG_GETARG_LINE_P(1); + + PG_RETURN_BOOL(lseg_interpt_line(NULL, lseg, line)); +} + + +/* + * Do line segment and box intersect? + * + * Segment completely inside box counts as intersection. + * If you want only segments crossing box boundaries, + * try converting box to path first. + * + * This function also sets the *result to the closest point on the line + * segment to the center of the box when they overlap and the result is + * not NULL. It is somewhat arbitrary, but maybe the best we can do as + * there are typically two points they intersect. + * + * Optimize for non-intersection by checking for box intersection first. + * - thomas 1998-01-30 + */ +static bool +box_interpt_lseg(Point *result, BOX *box, LSEG *lseg) +{ + BOX lbox; + LSEG bseg; + Point point; + + lbox.low.x = float8_min(lseg->p[0].x, lseg->p[1].x); + lbox.low.y = float8_min(lseg->p[0].y, lseg->p[1].y); + lbox.high.x = float8_max(lseg->p[0].x, lseg->p[1].x); + lbox.high.y = float8_max(lseg->p[0].y, lseg->p[1].y); + + /* nothing close to overlap? then not going to intersect */ + if (!box_ov(&lbox, box)) + return false; + + if (result != NULL) + { + box_cn(&point, box); + lseg_closept_point(result, lseg, &point); + } + + /* an endpoint of segment is inside box? then clearly intersects */ + if (box_contain_point(box, &lseg->p[0]) || + box_contain_point(box, &lseg->p[1])) + return true; + + /* pairwise check lseg intersections */ + point.x = box->low.x; + point.y = box->high.y; + statlseg_construct(&bseg, &box->low, &point); + if (lseg_interpt_lseg(NULL, &bseg, lseg)) + return true; + + statlseg_construct(&bseg, &box->high, &point); + if (lseg_interpt_lseg(NULL, &bseg, lseg)) + return true; + + point.x = box->high.x; + point.y = box->low.y; + statlseg_construct(&bseg, &box->low, &point); + if (lseg_interpt_lseg(NULL, &bseg, lseg)) + return true; + + statlseg_construct(&bseg, &box->high, &point); + if (lseg_interpt_lseg(NULL, &bseg, lseg)) + return true; + + /* if we dropped through, no two segs intersected */ + return false; +} + +Datum +inter_sb(PG_FUNCTION_ARGS) +{ + LSEG *lseg = PG_GETARG_LSEG_P(0); + BOX *box = PG_GETARG_BOX_P(1); + + PG_RETURN_BOOL(box_interpt_lseg(NULL, box, lseg)); +} + + +/* inter_lb() + * Do line and box intersect? + */ +Datum +inter_lb(PG_FUNCTION_ARGS) +{ + LINE *line = PG_GETARG_LINE_P(0); + BOX *box = PG_GETARG_BOX_P(1); + LSEG bseg; + Point p1, + p2; + + /* pairwise check lseg intersections */ + p1.x = box->low.x; + p1.y = box->low.y; + p2.x = box->low.x; + p2.y = box->high.y; + statlseg_construct(&bseg, &p1, &p2); + if (lseg_interpt_line(NULL, &bseg, line)) + PG_RETURN_BOOL(true); + p1.x = box->high.x; + p1.y = box->high.y; + statlseg_construct(&bseg, &p1, &p2); + if (lseg_interpt_line(NULL, &bseg, line)) + PG_RETURN_BOOL(true); + p2.x = box->high.x; + p2.y = box->low.y; + statlseg_construct(&bseg, &p1, &p2); + if (lseg_interpt_line(NULL, &bseg, line)) + PG_RETURN_BOOL(true); + p1.x = box->low.x; + p1.y = box->low.y; + statlseg_construct(&bseg, &p1, &p2); + if (lseg_interpt_line(NULL, &bseg, line)) + PG_RETURN_BOOL(true); + + /* if we dropped through, no intersection */ + PG_RETURN_BOOL(false); +} + +/*------------------------------------------------------------------ + * The following routines define a data type and operator class for + * POLYGONS .... Part of which (the polygon's bounding box) is built on + * top of the BOX data type. + * + * make_bound_box - create the bounding box for the input polygon + *------------------------------------------------------------------*/ + +/*--------------------------------------------------------------------- + * Make the smallest bounding box for the given polygon. + *---------------------------------------------------------------------*/ +static void +make_bound_box(POLYGON *poly) +{ + int i; + float8 x1, + y1, + x2, + y2; + + Assert(poly->npts > 0); + + x1 = x2 = poly->p[0].x; + y2 = y1 = poly->p[0].y; + for (i = 1; i < poly->npts; i++) + { + if (float8_lt(poly->p[i].x, x1)) + x1 = poly->p[i].x; + if (float8_gt(poly->p[i].x, x2)) + x2 = poly->p[i].x; + if (float8_lt(poly->p[i].y, y1)) + y1 = poly->p[i].y; + if (float8_gt(poly->p[i].y, y2)) + y2 = poly->p[i].y; + } + + poly->boundbox.low.x = x1; + poly->boundbox.high.x = x2; + poly->boundbox.low.y = y1; + poly->boundbox.high.y = y2; +} + +/*------------------------------------------------------------------ + * poly_in - read in the polygon from a string specification + * + * External format: + * "((x0,y0),...,(xn,yn))" + * "x0,y0,...,xn,yn" + * also supports the older style "(x1,...,xn,y1,...yn)" + *------------------------------------------------------------------*/ +Datum +poly_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + POLYGON *poly; + int npts; + int size; + int base_size; + bool isopen; + + if ((npts = pair_count(str, ',')) <= 0) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + "polygon", str))); + + base_size = sizeof(poly->p[0]) * npts; + size = offsetof(POLYGON, p) + base_size; + + /* Check for integer overflow */ + if (base_size / npts != sizeof(poly->p[0]) || size <= base_size) + ereport(ERROR, + (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED), + errmsg("too many points requested"))); + + poly = (POLYGON *) palloc0(size); /* zero any holes */ + + SET_VARSIZE(poly, size); + poly->npts = npts; + + path_decode(str, false, npts, &(poly->p[0]), &isopen, NULL, "polygon", str); + + make_bound_box(poly); + + PG_RETURN_POLYGON_P(poly); +} + +/*--------------------------------------------------------------- + * poly_out - convert internal POLYGON representation to the + * character string format "((f8,f8),...,(f8,f8))" + *---------------------------------------------------------------*/ +Datum +poly_out(PG_FUNCTION_ARGS) +{ + POLYGON *poly = PG_GETARG_POLYGON_P(0); + + PG_RETURN_CSTRING(path_encode(PATH_CLOSED, poly->npts, poly->p)); +} + +/* + * poly_recv - converts external binary format to polygon + * + * External representation is int32 number of points, and the points. + * We recompute the bounding box on read, instead of trusting it to + * be valid. (Checking it would take just as long, so may as well + * omit it from external representation.) + */ +Datum +poly_recv(PG_FUNCTION_ARGS) +{ + StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); + POLYGON *poly; + int32 npts; + int32 i; + int size; + + npts = pq_getmsgint(buf, sizeof(int32)); + if (npts <= 0 || npts >= (int32) ((INT_MAX - offsetof(POLYGON, p)) / sizeof(Point))) + ereport(ERROR, + (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION), + errmsg("invalid number of points in external \"polygon\" value"))); + + size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * npts; + poly = (POLYGON *) palloc0(size); /* zero any holes */ + + SET_VARSIZE(poly, size); + poly->npts = npts; + + for (i = 0; i < npts; i++) + { + poly->p[i].x = pq_getmsgfloat8(buf); + poly->p[i].y = pq_getmsgfloat8(buf); + } + + make_bound_box(poly); + + PG_RETURN_POLYGON_P(poly); +} + +/* + * poly_send - converts polygon to binary format + */ +Datum +poly_send(PG_FUNCTION_ARGS) +{ + POLYGON *poly = PG_GETARG_POLYGON_P(0); + StringInfoData buf; + int32 i; + + pq_begintypsend(&buf); + pq_sendint32(&buf, poly->npts); + for (i = 0; i < poly->npts; i++) + { + pq_sendfloat8(&buf, poly->p[i].x); + pq_sendfloat8(&buf, poly->p[i].y); + } + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); +} + + +/*------------------------------------------------------- + * Is polygon A strictly left of polygon B? i.e. is + * the right most point of A left of the left most point + * of B? + *-------------------------------------------------------*/ +Datum +poly_left(PG_FUNCTION_ARGS) +{ + POLYGON *polya = PG_GETARG_POLYGON_P(0); + POLYGON *polyb = PG_GETARG_POLYGON_P(1); + bool result; + + result = polya->boundbox.high.x < polyb->boundbox.low.x; + + /* + * Avoid leaking memory for toasted inputs ... needed for rtree indexes + */ + PG_FREE_IF_COPY(polya, 0); + PG_FREE_IF_COPY(polyb, 1); + + PG_RETURN_BOOL(result); +} + +/*------------------------------------------------------- + * Is polygon A overlapping or left of polygon B? i.e. is + * the right most point of A at or left of the right most point + * of B? + *-------------------------------------------------------*/ +Datum +poly_overleft(PG_FUNCTION_ARGS) +{ + POLYGON *polya = PG_GETARG_POLYGON_P(0); + POLYGON *polyb = PG_GETARG_POLYGON_P(1); + bool result; + + result = polya->boundbox.high.x <= polyb->boundbox.high.x; + + /* + * Avoid leaking memory for toasted inputs ... needed for rtree indexes + */ + PG_FREE_IF_COPY(polya, 0); + PG_FREE_IF_COPY(polyb, 1); + + PG_RETURN_BOOL(result); +} + +/*------------------------------------------------------- + * Is polygon A strictly right of polygon B? i.e. is + * the left most point of A right of the right most point + * of B? + *-------------------------------------------------------*/ +Datum +poly_right(PG_FUNCTION_ARGS) +{ + POLYGON *polya = PG_GETARG_POLYGON_P(0); + POLYGON *polyb = PG_GETARG_POLYGON_P(1); + bool result; + + result = polya->boundbox.low.x > polyb->boundbox.high.x; + + /* + * Avoid leaking memory for toasted inputs ... needed for rtree indexes + */ + PG_FREE_IF_COPY(polya, 0); + PG_FREE_IF_COPY(polyb, 1); + + PG_RETURN_BOOL(result); +} + +/*------------------------------------------------------- + * Is polygon A overlapping or right of polygon B? i.e. is + * the left most point of A at or right of the left most point + * of B? + *-------------------------------------------------------*/ +Datum +poly_overright(PG_FUNCTION_ARGS) +{ + POLYGON *polya = PG_GETARG_POLYGON_P(0); + POLYGON *polyb = PG_GETARG_POLYGON_P(1); + bool result; + + result = polya->boundbox.low.x >= polyb->boundbox.low.x; + + /* + * Avoid leaking memory for toasted inputs ... needed for rtree indexes + */ + PG_FREE_IF_COPY(polya, 0); + PG_FREE_IF_COPY(polyb, 1); + + PG_RETURN_BOOL(result); +} + +/*------------------------------------------------------- + * Is polygon A strictly below polygon B? i.e. is + * the upper most point of A below the lower most point + * of B? + *-------------------------------------------------------*/ +Datum +poly_below(PG_FUNCTION_ARGS) +{ + POLYGON *polya = PG_GETARG_POLYGON_P(0); + POLYGON *polyb = PG_GETARG_POLYGON_P(1); + bool result; + + result = polya->boundbox.high.y < polyb->boundbox.low.y; + + /* + * Avoid leaking memory for toasted inputs ... needed for rtree indexes + */ + PG_FREE_IF_COPY(polya, 0); + PG_FREE_IF_COPY(polyb, 1); + + PG_RETURN_BOOL(result); +} + +/*------------------------------------------------------- + * Is polygon A overlapping or below polygon B? i.e. is + * the upper most point of A at or below the upper most point + * of B? + *-------------------------------------------------------*/ +Datum +poly_overbelow(PG_FUNCTION_ARGS) +{ + POLYGON *polya = PG_GETARG_POLYGON_P(0); + POLYGON *polyb = PG_GETARG_POLYGON_P(1); + bool result; + + result = polya->boundbox.high.y <= polyb->boundbox.high.y; + + /* + * Avoid leaking memory for toasted inputs ... needed for rtree indexes + */ + PG_FREE_IF_COPY(polya, 0); + PG_FREE_IF_COPY(polyb, 1); + + PG_RETURN_BOOL(result); +} + +/*------------------------------------------------------- + * Is polygon A strictly above polygon B? i.e. is + * the lower most point of A above the upper most point + * of B? + *-------------------------------------------------------*/ +Datum +poly_above(PG_FUNCTION_ARGS) +{ + POLYGON *polya = PG_GETARG_POLYGON_P(0); + POLYGON *polyb = PG_GETARG_POLYGON_P(1); + bool result; + + result = polya->boundbox.low.y > polyb->boundbox.high.y; + + /* + * Avoid leaking memory for toasted inputs ... needed for rtree indexes + */ + PG_FREE_IF_COPY(polya, 0); + PG_FREE_IF_COPY(polyb, 1); + + PG_RETURN_BOOL(result); +} + +/*------------------------------------------------------- + * Is polygon A overlapping or above polygon B? i.e. is + * the lower most point of A at or above the lower most point + * of B? + *-------------------------------------------------------*/ +Datum +poly_overabove(PG_FUNCTION_ARGS) +{ + POLYGON *polya = PG_GETARG_POLYGON_P(0); + POLYGON *polyb = PG_GETARG_POLYGON_P(1); + bool result; + + result = polya->boundbox.low.y >= polyb->boundbox.low.y; + + /* + * Avoid leaking memory for toasted inputs ... needed for rtree indexes + */ + PG_FREE_IF_COPY(polya, 0); + PG_FREE_IF_COPY(polyb, 1); + + PG_RETURN_BOOL(result); +} + + +/*------------------------------------------------------- + * Is polygon A the same as polygon B? i.e. are all the + * points the same? + * Check all points for matches in both forward and reverse + * direction since polygons are non-directional and are + * closed shapes. + *-------------------------------------------------------*/ +Datum +poly_same(PG_FUNCTION_ARGS) +{ + POLYGON *polya = PG_GETARG_POLYGON_P(0); + POLYGON *polyb = PG_GETARG_POLYGON_P(1); + bool result; + + if (polya->npts != polyb->npts) + result = false; + else + result = plist_same(polya->npts, polya->p, polyb->p); + + /* + * Avoid leaking memory for toasted inputs ... needed for rtree indexes + */ + PG_FREE_IF_COPY(polya, 0); + PG_FREE_IF_COPY(polyb, 1); + + PG_RETURN_BOOL(result); +} + +/*----------------------------------------------------------------- + * Determine if polygon A overlaps polygon B + *-----------------------------------------------------------------*/ +Datum +poly_overlap(PG_FUNCTION_ARGS) +{ + POLYGON *polya = PG_GETARG_POLYGON_P(0); + POLYGON *polyb = PG_GETARG_POLYGON_P(1); + bool result; + + Assert(polya->npts > 0 && polyb->npts > 0); + + /* Quick check by bounding box */ + result = box_ov(&polya->boundbox, &polyb->boundbox); + + /* + * Brute-force algorithm - try to find intersected edges, if so then + * polygons are overlapped else check is one polygon inside other or not + * by testing single point of them. + */ + if (result) + { + int ia, + ib; + LSEG sa, + sb; + + /* Init first of polya's edge with last point */ + sa.p[0] = polya->p[polya->npts - 1]; + result = false; + + for (ia = 0; ia < polya->npts && !result; ia++) + { + /* Second point of polya's edge is a current one */ + sa.p[1] = polya->p[ia]; + + /* Init first of polyb's edge with last point */ + sb.p[0] = polyb->p[polyb->npts - 1]; + + for (ib = 0; ib < polyb->npts && !result; ib++) + { + sb.p[1] = polyb->p[ib]; + result = lseg_interpt_lseg(NULL, &sa, &sb); + sb.p[0] = sb.p[1]; + } + + /* + * move current endpoint to the first point of next edge + */ + sa.p[0] = sa.p[1]; + } + + if (!result) + { + result = (point_inside(polya->p, polyb->npts, polyb->p) || + point_inside(polyb->p, polya->npts, polya->p)); + } + } + + /* + * Avoid leaking memory for toasted inputs ... needed for rtree indexes + */ + PG_FREE_IF_COPY(polya, 0); + PG_FREE_IF_COPY(polyb, 1); + + PG_RETURN_BOOL(result); +} + +/* + * Tests special kind of segment for in/out of polygon. + * Special kind means: + * - point a should be on segment s + * - segment (a,b) should not be contained by s + * Returns true if: + * - segment (a,b) is collinear to s and (a,b) is in polygon + * - segment (a,b) s not collinear to s. Note: that doesn't + * mean that segment is in polygon! + */ + +static bool +touched_lseg_inside_poly(Point *a, Point *b, LSEG *s, POLYGON *poly, int start) +{ + /* point a is on s, b is not */ + LSEG t; + + t.p[0] = *a; + t.p[1] = *b; + + if (point_eq_point(a, s->p)) + { + if (lseg_contain_point(&t, s->p + 1)) + return lseg_inside_poly(b, s->p + 1, poly, start); + } + else if (point_eq_point(a, s->p + 1)) + { + if (lseg_contain_point(&t, s->p)) + return lseg_inside_poly(b, s->p, poly, start); + } + else if (lseg_contain_point(&t, s->p)) + { + return lseg_inside_poly(b, s->p, poly, start); + } + else if (lseg_contain_point(&t, s->p + 1)) + { + return lseg_inside_poly(b, s->p + 1, poly, start); + } + + return true; /* may be not true, but that will check later */ +} + +/* + * Returns true if segment (a,b) is in polygon, option + * start is used for optimization - function checks + * polygon's edges starting from start + */ +static bool +lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start) +{ + LSEG s, + t; + int i; + bool res = true, + intersection = false; + + t.p[0] = *a; + t.p[1] = *b; + s.p[0] = poly->p[(start == 0) ? (poly->npts - 1) : (start - 1)]; + + for (i = start; i < poly->npts && res; i++) + { + Point interpt; + + CHECK_FOR_INTERRUPTS(); + + s.p[1] = poly->p[i]; + + if (lseg_contain_point(&s, t.p)) + { + if (lseg_contain_point(&s, t.p + 1)) + return true; /* t is contained by s */ + + /* Y-cross */ + res = touched_lseg_inside_poly(t.p, t.p + 1, &s, poly, i + 1); + } + else if (lseg_contain_point(&s, t.p + 1)) + { + /* Y-cross */ + res = touched_lseg_inside_poly(t.p + 1, t.p, &s, poly, i + 1); + } + else if (lseg_interpt_lseg(&interpt, &t, &s)) + { + /* + * segments are X-crossing, go to check each subsegment + */ + + intersection = true; + res = lseg_inside_poly(t.p, &interpt, poly, i + 1); + if (res) + res = lseg_inside_poly(t.p + 1, &interpt, poly, i + 1); + } + + s.p[0] = s.p[1]; + } + + if (res && !intersection) + { + Point p; + + /* + * if X-intersection wasn't found then check central point of tested + * segment. In opposite case we already check all subsegments + */ + p.x = float8_div(float8_pl(t.p[0].x, t.p[1].x), 2.0); + p.y = float8_div(float8_pl(t.p[0].y, t.p[1].y), 2.0); + + res = point_inside(&p, poly->npts, poly->p); + } + + return res; +} + +/* + * Check whether the first polygon contains the second + */ +static bool +poly_contain_poly(POLYGON *contains_poly, POLYGON *contained_poly) +{ + int i; + LSEG s; + + Assert(contains_poly->npts > 0 && contained_poly->npts > 0); + + /* + * Quick check to see if contained's bounding box is contained in + * contains' bb. + */ + if (!box_contain_box(&contains_poly->boundbox, &contained_poly->boundbox)) + return false; + + s.p[0] = contained_poly->p[contained_poly->npts - 1]; + + for (i = 0; i < contained_poly->npts; i++) + { + s.p[1] = contained_poly->p[i]; + if (!lseg_inside_poly(s.p, s.p + 1, contains_poly, 0)) + return false; + s.p[0] = s.p[1]; + } + + return true; +} + +Datum +poly_contain(PG_FUNCTION_ARGS) +{ + POLYGON *polya = PG_GETARG_POLYGON_P(0); + POLYGON *polyb = PG_GETARG_POLYGON_P(1); + bool result; + + result = poly_contain_poly(polya, polyb); + + /* + * Avoid leaking memory for toasted inputs ... needed for rtree indexes + */ + PG_FREE_IF_COPY(polya, 0); + PG_FREE_IF_COPY(polyb, 1); + + PG_RETURN_BOOL(result); +} + + +/*----------------------------------------------------------------- + * Determine if polygon A is contained by polygon B + *-----------------------------------------------------------------*/ +Datum +poly_contained(PG_FUNCTION_ARGS) +{ + POLYGON *polya = PG_GETARG_POLYGON_P(0); + POLYGON *polyb = PG_GETARG_POLYGON_P(1); + bool result; + + /* Just switch the arguments and pass it off to poly_contain */ + result = poly_contain_poly(polyb, polya); + + /* + * Avoid leaking memory for toasted inputs ... needed for rtree indexes + */ + PG_FREE_IF_COPY(polya, 0); + PG_FREE_IF_COPY(polyb, 1); + + PG_RETURN_BOOL(result); +} + + +Datum +poly_contain_pt(PG_FUNCTION_ARGS) +{ + POLYGON *poly = PG_GETARG_POLYGON_P(0); + Point *p = PG_GETARG_POINT_P(1); + + PG_RETURN_BOOL(point_inside(p, poly->npts, poly->p) != 0); +} + +Datum +pt_contained_poly(PG_FUNCTION_ARGS) +{ + Point *p = PG_GETARG_POINT_P(0); + POLYGON *poly = PG_GETARG_POLYGON_P(1); + + PG_RETURN_BOOL(point_inside(p, poly->npts, poly->p) != 0); +} + + +Datum +poly_distance(PG_FUNCTION_ARGS) +{ +#ifdef NOT_USED + POLYGON *polya = PG_GETARG_POLYGON_P(0); + POLYGON *polyb = PG_GETARG_POLYGON_P(1); +#endif + + ereport(ERROR, + (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), + errmsg("function \"poly_distance\" not implemented"))); + + PG_RETURN_NULL(); +} + + +/*********************************************************************** + ** + ** Routines for 2D points. + ** + ***********************************************************************/ + +Datum +construct_point(PG_FUNCTION_ARGS) +{ + float8 x = PG_GETARG_FLOAT8(0); + float8 y = PG_GETARG_FLOAT8(1); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + point_construct(result, x, y); + + PG_RETURN_POINT_P(result); +} + + +static inline void +point_add_point(Point *result, Point *pt1, Point *pt2) +{ + point_construct(result, + float8_pl(pt1->x, pt2->x), + float8_pl(pt1->y, pt2->y)); +} + +Datum +point_add(PG_FUNCTION_ARGS) +{ + Point *p1 = PG_GETARG_POINT_P(0); + Point *p2 = PG_GETARG_POINT_P(1); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + point_add_point(result, p1, p2); + + PG_RETURN_POINT_P(result); +} + + +static inline void +point_sub_point(Point *result, Point *pt1, Point *pt2) +{ + point_construct(result, + float8_mi(pt1->x, pt2->x), + float8_mi(pt1->y, pt2->y)); +} + +Datum +point_sub(PG_FUNCTION_ARGS) +{ + Point *p1 = PG_GETARG_POINT_P(0); + Point *p2 = PG_GETARG_POINT_P(1); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + point_sub_point(result, p1, p2); + + PG_RETURN_POINT_P(result); +} + + +static inline void +point_mul_point(Point *result, Point *pt1, Point *pt2) +{ + point_construct(result, + float8_mi(float8_mul(pt1->x, pt2->x), + float8_mul(pt1->y, pt2->y)), + float8_pl(float8_mul(pt1->x, pt2->y), + float8_mul(pt1->y, pt2->x))); +} + +Datum +point_mul(PG_FUNCTION_ARGS) +{ + Point *p1 = PG_GETARG_POINT_P(0); + Point *p2 = PG_GETARG_POINT_P(1); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + point_mul_point(result, p1, p2); + + PG_RETURN_POINT_P(result); +} + + +static inline void +point_div_point(Point *result, Point *pt1, Point *pt2) +{ + float8 div; + + div = float8_pl(float8_mul(pt2->x, pt2->x), float8_mul(pt2->y, pt2->y)); + + point_construct(result, + float8_div(float8_pl(float8_mul(pt1->x, pt2->x), + float8_mul(pt1->y, pt2->y)), div), + float8_div(float8_mi(float8_mul(pt1->y, pt2->x), + float8_mul(pt1->x, pt2->y)), div)); +} + +Datum +point_div(PG_FUNCTION_ARGS) +{ + Point *p1 = PG_GETARG_POINT_P(0); + Point *p2 = PG_GETARG_POINT_P(1); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + point_div_point(result, p1, p2); + + PG_RETURN_POINT_P(result); +} + + +/*********************************************************************** + ** + ** Routines for 2D boxes. + ** + ***********************************************************************/ + +Datum +points_box(PG_FUNCTION_ARGS) +{ + Point *p1 = PG_GETARG_POINT_P(0); + Point *p2 = PG_GETARG_POINT_P(1); + BOX *result; + + result = (BOX *) palloc(sizeof(BOX)); + + box_construct(result, p1, p2); + + PG_RETURN_BOX_P(result); +} + +Datum +box_add(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + Point *p = PG_GETARG_POINT_P(1); + BOX *result; + + result = (BOX *) palloc(sizeof(BOX)); + + point_add_point(&result->high, &box->high, p); + point_add_point(&result->low, &box->low, p); + + PG_RETURN_BOX_P(result); +} + +Datum +box_sub(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + Point *p = PG_GETARG_POINT_P(1); + BOX *result; + + result = (BOX *) palloc(sizeof(BOX)); + + point_sub_point(&result->high, &box->high, p); + point_sub_point(&result->low, &box->low, p); + + PG_RETURN_BOX_P(result); +} + +Datum +box_mul(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + Point *p = PG_GETARG_POINT_P(1); + BOX *result; + Point high, + low; + + result = (BOX *) palloc(sizeof(BOX)); + + point_mul_point(&high, &box->high, p); + point_mul_point(&low, &box->low, p); + + box_construct(result, &high, &low); + + PG_RETURN_BOX_P(result); +} + +Datum +box_div(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + Point *p = PG_GETARG_POINT_P(1); + BOX *result; + Point high, + low; + + result = (BOX *) palloc(sizeof(BOX)); + + point_div_point(&high, &box->high, p); + point_div_point(&low, &box->low, p); + + box_construct(result, &high, &low); + + PG_RETURN_BOX_P(result); +} + +/* + * Convert point to empty box + */ +Datum +point_box(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + BOX *box; + + box = (BOX *) palloc(sizeof(BOX)); + + box->high.x = pt->x; + box->low.x = pt->x; + box->high.y = pt->y; + box->low.y = pt->y; + + PG_RETURN_BOX_P(box); +} + +/* + * Smallest bounding box that includes both of the given boxes + */ +Datum +boxes_bound_box(PG_FUNCTION_ARGS) +{ + BOX *box1 = PG_GETARG_BOX_P(0), + *box2 = PG_GETARG_BOX_P(1), + *container; + + container = (BOX *) palloc(sizeof(BOX)); + + container->high.x = float8_max(box1->high.x, box2->high.x); + container->low.x = float8_min(box1->low.x, box2->low.x); + container->high.y = float8_max(box1->high.y, box2->high.y); + container->low.y = float8_min(box1->low.y, box2->low.y); + + PG_RETURN_BOX_P(container); +} + + +/*********************************************************************** + ** + ** Routines for 2D paths. + ** + ***********************************************************************/ + +/* path_add() + * Concatenate two paths (only if they are both open). + */ +Datum +path_add(PG_FUNCTION_ARGS) +{ + PATH *p1 = PG_GETARG_PATH_P(0); + PATH *p2 = PG_GETARG_PATH_P(1); + PATH *result; + int size, + base_size; + int i; + + if (p1->closed || p2->closed) + PG_RETURN_NULL(); + + base_size = sizeof(p1->p[0]) * (p1->npts + p2->npts); + size = offsetof(PATH, p) + base_size; + + /* Check for integer overflow */ + if (base_size / sizeof(p1->p[0]) != (p1->npts + p2->npts) || + size <= base_size) + ereport(ERROR, + (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED), + errmsg("too many points requested"))); + + result = (PATH *) palloc(size); + + SET_VARSIZE(result, size); + result->npts = (p1->npts + p2->npts); + result->closed = p1->closed; + /* prevent instability in unused pad bytes */ + result->dummy = 0; + + for (i = 0; i < p1->npts; i++) + { + result->p[i].x = p1->p[i].x; + result->p[i].y = p1->p[i].y; + } + for (i = 0; i < p2->npts; i++) + { + result->p[i + p1->npts].x = p2->p[i].x; + result->p[i + p1->npts].y = p2->p[i].y; + } + + PG_RETURN_PATH_P(result); +} + +/* path_add_pt() + * Translation operators. + */ +Datum +path_add_pt(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P_COPY(0); + Point *point = PG_GETARG_POINT_P(1); + int i; + + for (i = 0; i < path->npts; i++) + point_add_point(&path->p[i], &path->p[i], point); + + PG_RETURN_PATH_P(path); +} + +Datum +path_sub_pt(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P_COPY(0); + Point *point = PG_GETARG_POINT_P(1); + int i; + + for (i = 0; i < path->npts; i++) + point_sub_point(&path->p[i], &path->p[i], point); + + PG_RETURN_PATH_P(path); +} + +/* path_mul_pt() + * Rotation and scaling operators. + */ +Datum +path_mul_pt(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P_COPY(0); + Point *point = PG_GETARG_POINT_P(1); + int i; + + for (i = 0; i < path->npts; i++) + point_mul_point(&path->p[i], &path->p[i], point); + + PG_RETURN_PATH_P(path); +} + +Datum +path_div_pt(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P_COPY(0); + Point *point = PG_GETARG_POINT_P(1); + int i; + + for (i = 0; i < path->npts; i++) + point_div_point(&path->p[i], &path->p[i], point); + + PG_RETURN_PATH_P(path); +} + + +Datum +path_center(PG_FUNCTION_ARGS) +{ +#ifdef NOT_USED + PATH *path = PG_GETARG_PATH_P(0); +#endif + + ereport(ERROR, + (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), + errmsg("function \"path_center\" not implemented"))); + + PG_RETURN_NULL(); +} + +Datum +path_poly(PG_FUNCTION_ARGS) +{ + PATH *path = PG_GETARG_PATH_P(0); + POLYGON *poly; + int size; + int i; + + /* This is not very consistent --- other similar cases return NULL ... */ + if (!path->closed) + ereport(ERROR, + (errcode(ERRCODE_INVALID_PARAMETER_VALUE), + errmsg("open path cannot be converted to polygon"))); + + /* + * Never overflows: the old size fit in MaxAllocSize, and the new size is + * just a small constant larger. + */ + size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * path->npts; + poly = (POLYGON *) palloc(size); + + SET_VARSIZE(poly, size); + poly->npts = path->npts; + + for (i = 0; i < path->npts; i++) + { + poly->p[i].x = path->p[i].x; + poly->p[i].y = path->p[i].y; + } + + make_bound_box(poly); + + PG_RETURN_POLYGON_P(poly); +} + + +/*********************************************************************** + ** + ** Routines for 2D polygons. + ** + ***********************************************************************/ + +Datum +poly_npoints(PG_FUNCTION_ARGS) +{ + POLYGON *poly = PG_GETARG_POLYGON_P(0); + + PG_RETURN_INT32(poly->npts); +} + + +Datum +poly_center(PG_FUNCTION_ARGS) +{ + POLYGON *poly = PG_GETARG_POLYGON_P(0); + Point *result; + CIRCLE circle; + + result = (Point *) palloc(sizeof(Point)); + + poly_to_circle(&circle, poly); + *result = circle.center; + + PG_RETURN_POINT_P(result); +} + + +Datum +poly_box(PG_FUNCTION_ARGS) +{ + POLYGON *poly = PG_GETARG_POLYGON_P(0); + BOX *box; + + box = (BOX *) palloc(sizeof(BOX)); + *box = poly->boundbox; + + PG_RETURN_BOX_P(box); +} + + +/* box_poly() + * Convert a box to a polygon. + */ +Datum +box_poly(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + POLYGON *poly; + int size; + + /* map four corners of the box to a polygon */ + size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * 4; + poly = (POLYGON *) palloc(size); + + SET_VARSIZE(poly, size); + poly->npts = 4; + + poly->p[0].x = box->low.x; + poly->p[0].y = box->low.y; + poly->p[1].x = box->low.x; + poly->p[1].y = box->high.y; + poly->p[2].x = box->high.x; + poly->p[2].y = box->high.y; + poly->p[3].x = box->high.x; + poly->p[3].y = box->low.y; + + box_construct(&poly->boundbox, &box->high, &box->low); + + PG_RETURN_POLYGON_P(poly); +} + + +Datum +poly_path(PG_FUNCTION_ARGS) +{ + POLYGON *poly = PG_GETARG_POLYGON_P(0); + PATH *path; + int size; + int i; + + /* + * Never overflows: the old size fit in MaxAllocSize, and the new size is + * smaller by a small constant. + */ + size = offsetof(PATH, p) + sizeof(path->p[0]) * poly->npts; + path = (PATH *) palloc(size); + + SET_VARSIZE(path, size); + path->npts = poly->npts; + path->closed = true; + /* prevent instability in unused pad bytes */ + path->dummy = 0; + + for (i = 0; i < poly->npts; i++) + { + path->p[i].x = poly->p[i].x; + path->p[i].y = poly->p[i].y; + } + + PG_RETURN_PATH_P(path); +} + + +/*********************************************************************** + ** + ** Routines for circles. + ** + ***********************************************************************/ + +/*---------------------------------------------------------- + * Formatting and conversion routines. + *---------------------------------------------------------*/ + +/* circle_in - convert a string to internal form. + * + * External format: (center and radius of circle) + * "<(f8,f8),f8>" + * also supports quick entry style "f8,f8,f8" + */ +Datum +circle_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + CIRCLE *circle = (CIRCLE *) palloc(sizeof(CIRCLE)); + char *s, + *cp; + int depth = 0; + + s = str; + while (isspace((unsigned char) *s)) + s++; + if (*s == LDELIM_C) + depth++, s++; + else if (*s == LDELIM) + { + /* If there are two left parens, consume the first one */ + cp = (s + 1); + while (isspace((unsigned char) *cp)) + cp++; + if (*cp == LDELIM) + depth++, s = cp; + } + + /* pair_decode will consume parens around the pair, if any */ + pair_decode(s, &circle->center.x, &circle->center.y, &s, "circle", str); + + if (*s == DELIM) + s++; + + circle->radius = single_decode(s, &s, "circle", str); + /* We have to accept NaN. */ + if (circle->radius < 0.0) + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + "circle", str))); + + while (depth > 0) + { + if ((*s == RDELIM) || ((*s == RDELIM_C) && (depth == 1))) + { + depth--; + s++; + while (isspace((unsigned char) *s)) + s++; + } + else + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + "circle", str))); + } + + if (*s != '\0') + ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + "circle", str))); + + PG_RETURN_CIRCLE_P(circle); +} + +/* circle_out - convert a circle to external form. + */ +Datum +circle_out(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + StringInfoData str; + + initStringInfo(&str); + + appendStringInfoChar(&str, LDELIM_C); + appendStringInfoChar(&str, LDELIM); + pair_encode(circle->center.x, circle->center.y, &str); + appendStringInfoChar(&str, RDELIM); + appendStringInfoChar(&str, DELIM); + single_encode(circle->radius, &str); + appendStringInfoChar(&str, RDELIM_C); + + PG_RETURN_CSTRING(str.data); +} + +/* + * circle_recv - converts external binary format to circle + */ +Datum +circle_recv(PG_FUNCTION_ARGS) +{ + StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); + CIRCLE *circle; + + circle = (CIRCLE *) palloc(sizeof(CIRCLE)); + + circle->center.x = pq_getmsgfloat8(buf); + circle->center.y = pq_getmsgfloat8(buf); + circle->radius = pq_getmsgfloat8(buf); + + /* We have to accept NaN. */ + if (circle->radius < 0.0) + ereport(ERROR, + (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION), + errmsg("invalid radius in external \"circle\" value"))); + + PG_RETURN_CIRCLE_P(circle); +} + +/* + * circle_send - converts circle to binary format + */ +Datum +circle_send(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + StringInfoData buf; + + pq_begintypsend(&buf); + pq_sendfloat8(&buf, circle->center.x); + pq_sendfloat8(&buf, circle->center.y); + pq_sendfloat8(&buf, circle->radius); + PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); +} + + +/*---------------------------------------------------------- + * Relational operators for CIRCLEs. + * <, >, <=, >=, and == are based on circle area. + *---------------------------------------------------------*/ + +/* circles identical? + * + * We consider NaNs values to be equal to each other to let those circles + * to be found. + */ +Datum +circle_same(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(((isnan(circle1->radius) && isnan(circle1->radius)) || + FPeq(circle1->radius, circle2->radius)) && + point_eq_point(&circle1->center, &circle2->center)); +} + +/* circle_overlap - does circle1 overlap circle2? + */ +Datum +circle_overlap(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center), + float8_pl(circle1->radius, circle2->radius))); +} + +/* circle_overleft - is the right edge of circle1 at or left of + * the right edge of circle2? + */ +Datum +circle_overleft(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPle(float8_pl(circle1->center.x, circle1->radius), + float8_pl(circle2->center.x, circle2->radius))); +} + +/* circle_left - is circle1 strictly left of circle2? + */ +Datum +circle_left(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPlt(float8_pl(circle1->center.x, circle1->radius), + float8_mi(circle2->center.x, circle2->radius))); +} + +/* circle_right - is circle1 strictly right of circle2? + */ +Datum +circle_right(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPgt(float8_mi(circle1->center.x, circle1->radius), + float8_pl(circle2->center.x, circle2->radius))); +} + +/* circle_overright - is the left edge of circle1 at or right of + * the left edge of circle2? + */ +Datum +circle_overright(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPge(float8_mi(circle1->center.x, circle1->radius), + float8_mi(circle2->center.x, circle2->radius))); +} + +/* circle_contained - is circle1 contained by circle2? + */ +Datum +circle_contained(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center), + float8_mi(circle2->radius, circle1->radius))); +} + +/* circle_contain - does circle1 contain circle2? + */ +Datum +circle_contain(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center), + float8_mi(circle1->radius, circle2->radius))); +} + + +/* circle_below - is circle1 strictly below circle2? + */ +Datum +circle_below(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPlt(float8_pl(circle1->center.y, circle1->radius), + float8_mi(circle2->center.y, circle2->radius))); +} + +/* circle_above - is circle1 strictly above circle2? + */ +Datum +circle_above(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPgt(float8_mi(circle1->center.y, circle1->radius), + float8_pl(circle2->center.y, circle2->radius))); +} + +/* circle_overbelow - is the upper edge of circle1 at or below + * the upper edge of circle2? + */ +Datum +circle_overbelow(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPle(float8_pl(circle1->center.y, circle1->radius), + float8_pl(circle2->center.y, circle2->radius))); +} + +/* circle_overabove - is the lower edge of circle1 at or above + * the lower edge of circle2? + */ +Datum +circle_overabove(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPge(float8_mi(circle1->center.y, circle1->radius), + float8_mi(circle2->center.y, circle2->radius))); +} + + +/* circle_relop - is area(circle1) relop area(circle2), within + * our accuracy constraint? + */ +Datum +circle_eq(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPeq(circle_ar(circle1), circle_ar(circle2))); +} + +Datum +circle_ne(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPne(circle_ar(circle1), circle_ar(circle2))); +} + +Datum +circle_lt(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPlt(circle_ar(circle1), circle_ar(circle2))); +} + +Datum +circle_gt(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPgt(circle_ar(circle1), circle_ar(circle2))); +} + +Datum +circle_le(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPle(circle_ar(circle1), circle_ar(circle2))); +} + +Datum +circle_ge(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + + PG_RETURN_BOOL(FPge(circle_ar(circle1), circle_ar(circle2))); +} + + +/*---------------------------------------------------------- + * "Arithmetic" operators on circles. + *---------------------------------------------------------*/ + +/* circle_add_pt() + * Translation operator. + */ +Datum +circle_add_pt(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + Point *point = PG_GETARG_POINT_P(1); + CIRCLE *result; + + result = (CIRCLE *) palloc(sizeof(CIRCLE)); + + point_add_point(&result->center, &circle->center, point); + result->radius = circle->radius; + + PG_RETURN_CIRCLE_P(result); +} + +Datum +circle_sub_pt(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + Point *point = PG_GETARG_POINT_P(1); + CIRCLE *result; + + result = (CIRCLE *) palloc(sizeof(CIRCLE)); + + point_sub_point(&result->center, &circle->center, point); + result->radius = circle->radius; + + PG_RETURN_CIRCLE_P(result); +} + + +/* circle_mul_pt() + * Rotation and scaling operators. + */ +Datum +circle_mul_pt(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + Point *point = PG_GETARG_POINT_P(1); + CIRCLE *result; + + result = (CIRCLE *) palloc(sizeof(CIRCLE)); + + point_mul_point(&result->center, &circle->center, point); + result->radius = float8_mul(circle->radius, HYPOT(point->x, point->y)); + + PG_RETURN_CIRCLE_P(result); +} + +Datum +circle_div_pt(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + Point *point = PG_GETARG_POINT_P(1); + CIRCLE *result; + + result = (CIRCLE *) palloc(sizeof(CIRCLE)); + + point_div_point(&result->center, &circle->center, point); + result->radius = float8_div(circle->radius, HYPOT(point->x, point->y)); + + PG_RETURN_CIRCLE_P(result); +} + + +/* circle_area - returns the area of the circle. + */ +Datum +circle_area(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + + PG_RETURN_FLOAT8(circle_ar(circle)); +} + + +/* circle_diameter - returns the diameter of the circle. + */ +Datum +circle_diameter(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + + PG_RETURN_FLOAT8(float8_mul(circle->radius, 2.0)); +} + + +/* circle_radius - returns the radius of the circle. + */ +Datum +circle_radius(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + + PG_RETURN_FLOAT8(circle->radius); +} + + +/* circle_distance - returns the distance between + * two circles. + */ +Datum +circle_distance(PG_FUNCTION_ARGS) +{ + CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); + CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); + float8 result; + + result = float8_mi(point_dt(&circle1->center, &circle2->center), + float8_pl(circle1->radius, circle2->radius)); + if (result < 0.0) + result = 0.0; + + PG_RETURN_FLOAT8(result); +} + + +Datum +circle_contain_pt(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + Point *point = PG_GETARG_POINT_P(1); + float8 d; + + d = point_dt(&circle->center, point); + PG_RETURN_BOOL(d <= circle->radius); +} + + +Datum +pt_contained_circle(PG_FUNCTION_ARGS) +{ + Point *point = PG_GETARG_POINT_P(0); + CIRCLE *circle = PG_GETARG_CIRCLE_P(1); + float8 d; + + d = point_dt(&circle->center, point); + PG_RETURN_BOOL(d <= circle->radius); +} + + +/* dist_pc - returns the distance between + * a point and a circle. + */ +Datum +dist_pc(PG_FUNCTION_ARGS) +{ + Point *point = PG_GETARG_POINT_P(0); + CIRCLE *circle = PG_GETARG_CIRCLE_P(1); + float8 result; + + result = float8_mi(point_dt(point, &circle->center), + circle->radius); + if (result < 0.0) + result = 0.0; + + PG_RETURN_FLOAT8(result); +} + +/* + * Distance from a circle to a point + */ +Datum +dist_cpoint(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + Point *point = PG_GETARG_POINT_P(1); + float8 result; + + result = float8_mi(point_dt(point, &circle->center), circle->radius); + if (result < 0.0) + result = 0.0; + + PG_RETURN_FLOAT8(result); +} + +/* circle_center - returns the center point of the circle. + */ +Datum +circle_center(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + result->x = circle->center.x; + result->y = circle->center.y; + + PG_RETURN_POINT_P(result); +} + + +/* circle_ar - returns the area of the circle. + */ +static float8 +circle_ar(CIRCLE *circle) +{ + return float8_mul(float8_mul(circle->radius, circle->radius), M_PI); +} + + +/*---------------------------------------------------------- + * Conversion operators. + *---------------------------------------------------------*/ + +Datum +cr_circle(PG_FUNCTION_ARGS) +{ + Point *center = PG_GETARG_POINT_P(0); + float8 radius = PG_GETARG_FLOAT8(1); + CIRCLE *result; + + result = (CIRCLE *) palloc(sizeof(CIRCLE)); + + result->center.x = center->x; + result->center.y = center->y; + result->radius = radius; + + PG_RETURN_CIRCLE_P(result); +} + +Datum +circle_box(PG_FUNCTION_ARGS) +{ + CIRCLE *circle = PG_GETARG_CIRCLE_P(0); + BOX *box; + float8 delta; + + box = (BOX *) palloc(sizeof(BOX)); + + delta = float8_div(circle->radius, sqrt(2.0)); + + box->high.x = float8_pl(circle->center.x, delta); + box->low.x = float8_mi(circle->center.x, delta); + box->high.y = float8_pl(circle->center.y, delta); + box->low.y = float8_mi(circle->center.y, delta); + + PG_RETURN_BOX_P(box); +} + +/* box_circle() + * Convert a box to a circle. + */ +Datum +box_circle(PG_FUNCTION_ARGS) +{ + BOX *box = PG_GETARG_BOX_P(0); + CIRCLE *circle; + + circle = (CIRCLE *) palloc(sizeof(CIRCLE)); + + circle->center.x = float8_div(float8_pl(box->high.x, box->low.x), 2.0); + circle->center.y = float8_div(float8_pl(box->high.y, box->low.y), 2.0); + + circle->radius = point_dt(&circle->center, &box->high); + + PG_RETURN_CIRCLE_P(circle); +} + + +Datum +circle_poly(PG_FUNCTION_ARGS) +{ + int32 npts = PG_GETARG_INT32(0); + CIRCLE *circle = PG_GETARG_CIRCLE_P(1); + POLYGON *poly; + int base_size, + size; + int i; + float8 angle; + float8 anglestep; + + if (FPzero(circle->radius)) + ereport(ERROR, + (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), + errmsg("cannot convert circle with radius zero to polygon"))); + + if (npts < 2) + ereport(ERROR, + (errcode(ERRCODE_INVALID_PARAMETER_VALUE), + errmsg("must request at least 2 points"))); + + base_size = sizeof(poly->p[0]) * npts; + size = offsetof(POLYGON, p) + base_size; + + /* Check for integer overflow */ + if (base_size / npts != sizeof(poly->p[0]) || size <= base_size) + ereport(ERROR, + (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED), + errmsg("too many points requested"))); + + poly = (POLYGON *) palloc0(size); /* zero any holes */ + SET_VARSIZE(poly, size); + poly->npts = npts; + + anglestep = float8_div(2.0 * M_PI, npts); + + for (i = 0; i < npts; i++) + { + angle = float8_mul(anglestep, i); + + poly->p[i].x = float8_mi(circle->center.x, + float8_mul(circle->radius, cos(angle))); + poly->p[i].y = float8_pl(circle->center.y, + float8_mul(circle->radius, sin(angle))); + } + + make_bound_box(poly); + + PG_RETURN_POLYGON_P(poly); +} + +/* + * Convert polygon to circle + * + * The result must be preallocated. + * + * XXX This algorithm should use weighted means of line segments + * rather than straight average values of points - tgl 97/01/21. + */ +static void +poly_to_circle(CIRCLE *result, POLYGON *poly) +{ + int i; + + Assert(poly->npts > 0); + + result->center.x = 0; + result->center.y = 0; + result->radius = 0; + + for (i = 0; i < poly->npts; i++) + point_add_point(&result->center, &result->center, &poly->p[i]); + result->center.x = float8_div(result->center.x, poly->npts); + result->center.y = float8_div(result->center.y, poly->npts); + + for (i = 0; i < poly->npts; i++) + result->radius = float8_pl(result->radius, + point_dt(&poly->p[i], &result->center)); + result->radius = float8_div(result->radius, poly->npts); +} + +Datum +poly_circle(PG_FUNCTION_ARGS) +{ + POLYGON *poly = PG_GETARG_POLYGON_P(0); + CIRCLE *result; + + result = (CIRCLE *) palloc(sizeof(CIRCLE)); + + poly_to_circle(result, poly); + + PG_RETURN_CIRCLE_P(result); +} + + +/*********************************************************************** + ** + ** Private routines for multiple types. + ** + ***********************************************************************/ + +/* + * Test to see if the point is inside the polygon, returns 1/0, or 2 if + * the point is on the polygon. + * Code adapted but not copied from integer-based routines in WN: A + * Server for the HTTP + * version 1.15.1, file wn/image.c + * http://hopf.math.northwestern.edu/index.html + * Description of algorithm: http://www.linuxjournal.com/article/2197 + * http://www.linuxjournal.com/article/2029 + */ + +#define POINT_ON_POLYGON INT_MAX + +static int +point_inside(Point *p, int npts, Point *plist) +{ + float8 x0, + y0; + float8 prev_x, + prev_y; + int i = 0; + float8 x, + y; + int cross, + total_cross = 0; + + Assert(npts > 0); + + /* compute first polygon point relative to single point */ + x0 = float8_mi(plist[0].x, p->x); + y0 = float8_mi(plist[0].y, p->y); + + prev_x = x0; + prev_y = y0; + /* loop over polygon points and aggregate total_cross */ + for (i = 1; i < npts; i++) + { + /* compute next polygon point relative to single point */ + x = float8_mi(plist[i].x, p->x); + y = float8_mi(plist[i].y, p->y); + + /* compute previous to current point crossing */ + if ((cross = lseg_crossing(x, y, prev_x, prev_y)) == POINT_ON_POLYGON) + return 2; + total_cross += cross; + + prev_x = x; + prev_y = y; + } + + /* now do the first point */ + if ((cross = lseg_crossing(x0, y0, prev_x, prev_y)) == POINT_ON_POLYGON) + return 2; + total_cross += cross; + + if (total_cross != 0) + return 1; + return 0; +} + + +/* lseg_crossing() + * Returns +/-2 if line segment crosses the positive X-axis in a +/- direction. + * Returns +/-1 if one point is on the positive X-axis. + * Returns 0 if both points are on the positive X-axis, or there is no crossing. + * Returns POINT_ON_POLYGON if the segment contains (0,0). + * Wow, that is one confusing API, but it is used above, and when summed, + * can tell is if a point is in a polygon. + */ + +static int +lseg_crossing(float8 x, float8 y, float8 prev_x, float8 prev_y) +{ + float8 z; + int y_sign; + + if (FPzero(y)) + { /* y == 0, on X axis */ + if (FPzero(x)) /* (x,y) is (0,0)? */ + return POINT_ON_POLYGON; + else if (FPgt(x, 0)) + { /* x > 0 */ + if (FPzero(prev_y)) /* y and prev_y are zero */ + /* prev_x > 0? */ + return FPgt(prev_x, 0.0) ? 0 : POINT_ON_POLYGON; + return FPlt(prev_y, 0.0) ? 1 : -1; + } + else + { /* x < 0, x not on positive X axis */ + if (FPzero(prev_y)) + /* prev_x < 0? */ + return FPlt(prev_x, 0.0) ? 0 : POINT_ON_POLYGON; + return 0; + } + } + else + { /* y != 0 */ + /* compute y crossing direction from previous point */ + y_sign = FPgt(y, 0.0) ? 1 : -1; + + if (FPzero(prev_y)) + /* previous point was on X axis, so new point is either off or on */ + return FPlt(prev_x, 0.0) ? 0 : y_sign; + else if ((y_sign < 0 && FPlt(prev_y, 0.0)) || + (y_sign > 0 && FPgt(prev_y, 0.0))) + /* both above or below X axis */ + return 0; /* same sign */ + else + { /* y and prev_y cross X-axis */ + if (FPge(x, 0.0) && FPgt(prev_x, 0.0)) + /* both non-negative so cross positive X-axis */ + return 2 * y_sign; + if (FPlt(x, 0.0) && FPle(prev_x, 0.0)) + /* both non-positive so do not cross positive X-axis */ + return 0; + + /* x and y cross axes, see URL above point_inside() */ + z = float8_mi(float8_mul(float8_mi(x, prev_x), y), + float8_mul(float8_mi(y, prev_y), x)); + if (FPzero(z)) + return POINT_ON_POLYGON; + if ((y_sign < 0 && FPlt(z, 0.0)) || + (y_sign > 0 && FPgt(z, 0.0))) + return 0; + return 2 * y_sign; + } + } +} + + +static bool +plist_same(int npts, Point *p1, Point *p2) +{ + int i, + ii, + j; + + /* find match for first point */ + for (i = 0; i < npts; i++) + { + if (point_eq_point(&p2[i], &p1[0])) + { + + /* match found? then look forward through remaining points */ + for (ii = 1, j = i + 1; ii < npts; ii++, j++) + { + if (j >= npts) + j = 0; + if (!point_eq_point(&p2[j], &p1[ii])) + break; + } + if (ii == npts) + return true; + + /* match not found forwards? then look backwards */ + for (ii = 1, j = i - 1; ii < npts; ii++, j--) + { + if (j < 0) + j = (npts - 1); + if (!point_eq_point(&p2[j], &p1[ii])) + break; + } + if (ii == npts) + return true; + } + } + + return false; +} + + +/*------------------------------------------------------------------------- + * Determine the hypotenuse. + * + * If required, x and y are swapped to make x the larger number. The + * traditional formula of x^2+y^2 is rearranged to factor x outside the + * sqrt. This allows computation of the hypotenuse for significantly + * larger values, and with a higher precision than when using the naive + * formula. In particular, this cannot overflow unless the final result + * would be out-of-range. + * + * sqrt( x^2 + y^2 ) = sqrt( x^2( 1 + y^2/x^2) ) + * = x * sqrt( 1 + y^2/x^2 ) + * = x * sqrt( 1 + y/x * y/x ) + * + * It is expected that this routine will eventually be replaced with the + * C99 hypot() function. + * + * This implementation conforms to IEEE Std 1003.1 and GLIBC, in that the + * case of hypot(inf,nan) results in INF, and not NAN. + *----------------------------------------------------------------------- + */ +float8 +pg_hypot(float8 x, float8 y) +{ + float8 yx, + result; + + /* Handle INF and NaN properly */ + if (isinf(x) || isinf(y)) + return get_float8_infinity(); + + if (isnan(x) || isnan(y)) + return get_float8_nan(); + + /* Else, drop any minus signs */ + x = fabs(x); + y = fabs(y); + + /* Swap x and y if needed to make x the larger one */ + if (x < y) + { + float8 temp = x; + + x = y; + y = temp; + } + + /* + * If y is zero, the hypotenuse is x. This test saves a few cycles in + * such cases, but more importantly it also protects against + * divide-by-zero errors, since now x >= y. + */ + if (y == 0.0) + return x; + + /* Determine the hypotenuse */ + yx = y / x; + result = x * sqrt(1.0 + (yx * yx)); + + if (unlikely(isinf(result))) + float_overflow_error(); + if (unlikely(result == 0.0)) + float_underflow_error(); + + return result; +} |