From 096a1e74ba77515610764ce362aa3d2ab87b7073 Mon Sep 17 00:00:00 2001 From: Emre Hasegeli Date: Sat, 27 May 2017 11:27:18 -0400 Subject: [PATCH 1/6] geo-funcs-v08 Refactor geometric functions and operators code Most of the geometric types were not using functions of each other's. I believe the reason behind this is simpler types line point and line being developed after more complicated ones. This patch reduces duplicate code and makes functions of different datatypes more compatible. The changes can be summarised as: * Eliminate SQL level function calls * Re-use more functions to implement others * Unify internal function names and signatures * Remove private functions from geo_decls.h * Replace not-possible checks with Assert() * Add comments to describe some functions * Remove some unreachable code * Define delimiter symbols of line datatype like the other ones * Remove debugging print()s from the refactored functions This commits aims to not cause any user visible changes, but it is almost impossible to stay completely bug-compatible with the old code. Discussion: https://www.postgresql.org/message-id/CAE2gYzxF7-5djV6-cEvqQu-fNsnt%3DEqbOURx7ZDg%2BVv6ZMTWbg%40mail.gmail.com --- src/backend/utils/adt/geo_ops.c | 1906 +++++++++++++++++------------------- src/backend/utils/adt/geo_spgist.c | 3 +- src/include/utils/geo_decls.h | 4 - src/test/regress/regress.c | 7 +- 4 files changed, 880 insertions(+), 1040 deletions(-) diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c index f57380a4df..4f8ac71a6b 100644 --- a/src/backend/utils/adt/geo_ops.c +++ b/src/backend/utils/adt/geo_ops.c @@ -31,75 +31,104 @@ /* * 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_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(double x, double y, double px, double py); -static BOX *box_construct(double x1, double x2, double y1, double y2); -static BOX *box_fill(BOX *result, double x1, double x2, double y1, double y2); +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 *l1, LSEG *l2); + +/* 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 double box_ar(BOX *box); static double box_ht(BOX *box); static double box_wd(BOX *box); +static bool box_contain_point(BOX *box, Point *point); +static bool box_contain_box(BOX *box1, BOX *box2); +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 double circle_ar(CIRCLE *circle); -static CIRCLE *circle_copy(CIRCLE *circle); -static LINE *line_construct_pm(Point *pt, double m); -static void line_construct_pts(LINE *line, Point *pt1, Point *pt2); -static bool lseg_intersect_internal(LSEG *l1, LSEG *l2); -static double lseg_dt(LSEG *l1, LSEG *l2); -static bool on_ps_internal(Point *pt, LSEG *lseg); + +/* 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 *polya, POLYGON *polyb); static bool plist_same(int npts, Point *p1, Point *p2); -static Point *point_construct(double x, double y); -static Point *point_copy(Point *pt); +static float8 dist_ppoly_internal(Point *pt, POLYGON *poly); + +/* Routines for encoding and decoding */ static double 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, double *x, double *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); -static void statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2); -static double box_ar(BOX *box); -static void box_cn(Point *center, BOX *box); -static Point *interpt_sl(LSEG *lseg, LINE *line); -static bool has_interpt_sl(LSEG *lseg, LINE *line); -static double dist_pl_internal(Point *pt, LINE *line); -static double dist_ps_internal(Point *pt, LSEG *lseg); -static Point *line_interpt_internal(LINE *l1, LINE *l2); -static bool lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start); -static Point *lseg_interpt_internal(LSEG *l1, LSEG *l2); -static double dist_ppoly_internal(Point *pt, POLYGON *poly); /* * 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). * @@ -229,26 +258,23 @@ path_decode(char *str, bool opentype, int npts, Point *p, } for (i = 0; i < npts; i++) { pair_decode(str, &(p->x), &(p->y), &str, type_name, orig_string); if (*str == DELIM) str++; p++; } - while (isspace((unsigned char) *str)) - str++; while (depth > 0) { - if ((*str == RDELIM) - || ((*str == RDELIM_EP) && (*isopen) && (depth == 1))) + 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\"", @@ -433,89 +459,61 @@ box_send(PG_FUNCTION_ARGS) 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 BOX * -box_construct(double x1, double x2, double y1, double y2) +static inline void +box_construct(BOX *result, Point *pt1, Point *pt2) { - BOX *result = (BOX *) palloc(sizeof(BOX)); - - return box_fill(result, x1, x2, y1, y2); -} - - -/* box_fill - fill in a given box struct - */ -static BOX * -box_fill(BOX *result, double x1, double x2, double y1, double y2) -{ - if (x1 > x2) + if (pt1->x > pt2->x) { - result->high.x = x1; - result->low.x = x2; + result->high.x = pt1->x; + result->low.x = pt2->x; } else { - result->high.x = x2; - result->low.x = x1; + result->high.x = pt2->x; + result->low.x = pt1->x; } - if (y1 > y2) + if (pt1->y > pt2->y) { - result->high.y = y1; - result->low.y = y2; + result->high.y = pt1->y; + result->low.y = pt2->y; } else { - result->high.y = y2; - result->low.y = y1; + result->high.y = pt2->y; + result->low.y = pt1->y; } - - return result; -} - - -/* box_copy - copy a box - */ -BOX * -box_copy(BOX *box) -{ - BOX *result = (BOX *) palloc(sizeof(BOX)); - - memcpy((char *) result, (char *) box, sizeof(BOX)); - - return result; } /*---------------------------------------------------------- * 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(FPeq(box1->high.x, box2->high.x) && - FPeq(box1->low.x, box2->low.x) && - FPeq(box1->high.y, box2->high.y) && - FPeq(box1->low.y, box2->low.y)); + 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); @@ -630,38 +628,44 @@ box_overabove(PG_FUNCTION_ARGS) } /* 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(FPle(box1->high.x, box2->high.x) && - FPge(box1->low.x, box2->low.x) && - FPle(box1->high.y, box2->high.y) && - FPge(box1->low.y, box2->low.y)); + 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(FPge(box1->high.x, box2->high.x) && - FPle(box1->low.x, box2->low.x) && - FPge(box1->high.y, box2->high.y) && - FPle(box1->low.y, box2->low.y)); + PG_RETURN_BOOL(box_contain_box(box1, box2)); +} + +/* + * Check whether the box is in the box or on its border + */ +static bool +box_contain_box(BOX *box1, BOX *box2) +{ + return FPge(box1->high.x, box2->high.x) && + FPle(box1->low.x, box2->low.x) && + FPge(box1->high.y, box2->high.y) && + FPle(box1->low.y, box2->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. @@ -750,51 +754,51 @@ box_area(PG_FUNCTION_ARGS) /* 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->high.x - box->low.x); + 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->high.y - box->low.y); + 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(HYPOT(a.x - b.x, a.y - b.y)); + 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)); @@ -898,76 +902,77 @@ 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++ != '}') + 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 == '{') + 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 (FPeq(lseg.p[0].x, lseg.p[1].x) && FPeq(lseg.p[0].y, lseg.p[1].y)) + 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_pts(line, &lseg.p[0], &lseg.p[1]); + 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("{%s,%s,%s}", astr, bstr, cstr)); + 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; @@ -996,128 +1001,78 @@ line_send(PG_FUNCTION_ARGS) 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 *---------------------------------------------------------*/ -/* line_construct_pm() - * point-slope +/* + * Fill already-allocated LINE struct from the point and the slope */ -static LINE * -line_construct_pm(Point *pt, double m) +static inline void +line_construct(LINE *result, Point *pt, float8 m) { - LINE *result = (LINE *) palloc(sizeof(LINE)); - if (m == DBL_MAX) { /* vertical - use "x = C" */ result->A = -1; result->B = 0; result->C = pt->x; } else { /* use "mx - y + yinter = 0" */ result->A = m; result->B = -1.0; result->C = pt->y - m * pt->x; } - - return result; -} - -/* - * Fill already-allocated LINE struct from two points on the line - */ -static void -line_construct_pts(LINE *line, Point *pt1, Point *pt2) -{ - if (FPeq(pt1->x, pt2->x)) - { /* vertical */ - /* use "x = C" */ - line->A = -1; - line->B = 0; - line->C = pt1->x; -#ifdef GEODEBUG - printf("line_construct_pts- line is vertical\n"); -#endif - } - else if (FPeq(pt1->y, pt2->y)) - { /* horizontal */ - /* use "y = C" */ - line->A = 0; - line->B = -1; - line->C = pt1->y; -#ifdef GEODEBUG - printf("line_construct_pts- line is horizontal\n"); -#endif - } - else - { - /* use "mx - y + yinter = 0" */ - line->A = (pt2->y - pt1->y) / (pt2->x - pt1->x); - line->B = -1.0; - line->C = pt1->y - line->A * pt1->x; - /* on some platforms, the preceding expression tends to produce -0 */ - if (line->C == 0.0) - line->C = 0.0; -#ifdef GEODEBUG - printf("line_construct_pts- line is neither vertical nor horizontal (diffs x=%.*g, y=%.*g\n", - DBL_DIG, (pt2->x - pt1->x), DBL_DIG, (pt2->y - pt1->y)); -#endif - } } /* 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)); - line_construct_pts(result, pt1, pt2); + 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(!DatumGetBool(DirectFunctionCall2(line_parallel, - LinePGetDatum(l1), - LinePGetDatum(l2)))); + 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); - if (FPzero(l1->B)) - PG_RETURN_BOOL(FPzero(l2->B)); - - PG_RETURN_BOOL(FPeq(l2->A, l1->A * (l2->B / l1->B))); + 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)); @@ -1162,105 +1117,113 @@ line_eq(PG_FUNCTION_ARGS) PG_RETURN_BOOL(FPeq(l1->A, k * l2->A) && FPeq(l1->B, k * l2->B) && FPeq(l1->C, k * l2->C)); } /*---------------------------------------------------------- * Line arithmetic routines. *---------------------------------------------------------*/ +/* + * Return inverse slope of the line + */ +static inline float8 +line_invsl(LINE *line) +{ + if (FPzero(line->A)) + return DBL_MAX; + if (FPzero(line->B)) + return 0.0; + return 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 result; - Point *tmp; + Point tmp; - if (!DatumGetBool(DirectFunctionCall2(line_parallel, - LinePGetDatum(l1), - LinePGetDatum(l2)))) + if (line_interpt_line(NULL, l1, l2)) PG_RETURN_FLOAT8(0.0); if (FPzero(l1->B)) /* vertical? */ PG_RETURN_FLOAT8(fabs(l1->C - l2->C)); - tmp = point_construct(0.0, l1->C); - result = dist_pl_internal(tmp, l2); + point_construct(&tmp, 0.0, l1->C); + result = line_closept_point(NULL, l2, &tmp); PG_RETURN_FLOAT8(result); } /* 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 = line_interpt_internal(l1, l2); + result = (Point *) palloc(sizeof(Point)); - if (result == NULL) + if (!line_interpt_line(result, l1, l2)) PG_RETURN_NULL(); PG_RETURN_POINT_P(result); } /* * Internal version of line_interpt * - * returns a NULL pointer if no intersection point + * This returns true if two lines intersect (they do, if they are not + * parallel), false if they do not. This also sets the intersection point + * to *result, if it is not NULL. + * + * 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. */ -static Point * -line_interpt_internal(LINE *l1, LINE *l2) +static bool +line_interpt_line(Point *result, LINE *l1, LINE *l2) { - Point *result; double x, y; - /* - * 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 (DatumGetBool(DirectFunctionCall2(line_parallel, - LinePGetDatum(l1), - LinePGetDatum(l2)))) - return NULL; - if (FPzero(l1->B)) /* l1 vertical? */ { + if (FPzero(l2->B)) /* l2 vertical? */ + return false; + x = l1->C; y = (l2->A * x + l2->C); } else if (FPzero(l2->B)) /* l2 vertical? */ { x = l2->C; y = (l1->A * x + l1->C); } else { + if (FPeq(l2->A, l1->A * (l2->B / l1->B))) + return false; + x = (l1->C - l2->C) / (l2->A - l1->A); y = (l1->A * x + l1->C); } - result = point_construct(x, y); + if (result != NULL) + point_construct(result, x, y); -#ifdef GEODEBUG - printf("line_interpt- lines are A=%.*g, B=%.*g, C=%.*g, A=%.*g, B=%.*g, C=%.*g\n", - DBL_DIG, l1->A, DBL_DIG, l1->B, DBL_DIG, l1->C, DBL_DIG, l2->A, DBL_DIG, l2->B, DBL_DIG, l2->C); - printf("line_interpt- lines intersect at (%.*g,%.*g)\n", DBL_DIG, x, DBL_DIG, y); -#endif - - return result; + 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 @@ -1555,22 +1518,21 @@ 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; - if (p1->npts <= 0 || p2->npts <= 0) - PG_RETURN_BOOL(false); + 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 = Max(p1->p[i].x, b1.high.x); b1.high.y = Max(p1->p[i].y, b1.high.y); b1.low.x = Min(p1->p[i].x, b1.low.x); b1.low.y = Min(p1->p[i].y, b1.low.y); } @@ -1608,21 +1570,21 @@ path_inter(PG_FUNCTION_ARGS) 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_intersect_internal(&seg1, &seg2)) + 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 @@ -1663,23 +1625,21 @@ path_distance(PG_FUNCTION_ARGS) 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 = DatumGetFloat8(DirectFunctionCall2(lseg_distance, - LsegPGetDatum(&seg1), - LsegPGetDatum(&seg2))); + tmp = lseg_closept_lseg(NULL, &seg1, &seg2); if (!have_min || tmp < min) { min = tmp; have_min = true; } } } if (!have_min) PG_RETURN_NULL(); @@ -1773,44 +1733,31 @@ 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)); } -static Point * -point_construct(double x, double y) +/* + * Initialize a point + * + * This function is over-simple, but it is, at least, useful when the new + * point is calculated using the existing one. + */ +static inline void +point_construct(Point *result, float8 x, float8 y) { - Point *result = (Point *) palloc(sizeof(Point)); - result->x = x; result->y = y; - return result; -} - - -static Point * -point_copy(Point *pt) -{ - Point *result; - - if (!PointerIsValid(pt)) - return NULL; - - result = (Point *) palloc(sizeof(Point)); - - result->x = pt->x; - result->y = pt->y; - return result; } /*---------------------------------------------------------- * 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). @@ -1869,71 +1816,97 @@ point_horiz(PG_FUNCTION_ARGS) 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(FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y)); + 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(FPne(pt1->x, pt2->x) || FPne(pt1->y, pt2->y)); + PG_RETURN_BOOL(!point_eq_point(pt1, pt2)); } +static inline bool +point_eq_point(Point *pt1, Point *pt2) +{ + 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(HYPOT(pt1->x - pt2->x, pt1->y - pt2->y)); + PG_RETURN_FLOAT8(point_dt(pt1, pt2)); } -double +static inline float8 point_dt(Point *pt1, Point *pt2) { -#ifdef GEODEBUG - printf("point_dt- segment (%f,%f),(%f,%f) length is %f\n", - pt1->x, pt1->y, pt2->x, pt2->y, HYPOT(pt1->x - pt2->x, pt1->y - pt2->y)); -#endif return HYPOT(pt1->x - pt2->x, 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)); } -double +/* + * Return slope of two points + * + * Note that this function returns DBL_MAX when the points are the same. + */ +static inline float8 point_sl(Point *pt1, Point *pt2) { - return (FPeq(pt1->x, pt2->x) - ? (double) DBL_MAX - : (pt1->y - pt2->y) / (pt1->x - pt2->x)); + if (FPeq(pt1->x, pt2->x)) + return DBL_MAX; + if (FPeq(pt1->y, pt2->y)) + return 0.0; + return (pt1->y - pt2->y) / (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 DBL_MAX; + return (pt1->x - pt2->x) / (pt2->y - pt1->y); } /*********************************************************************** ** ** Routines for 2D line segments. ** ***********************************************************************/ /*---------------------------------------------------------- @@ -1945,31 +1918,31 @@ point_sl(Point *pt1, Point *pt2) * (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); + 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, (Point *) &(ls->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; @@ -2005,38 +1978,56 @@ lseg_send(PG_FUNCTION_ARGS) /* 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)); - result->p[0].x = pt1->x; - result->p[0].y = pt1->y; - result->p[1].x = pt2->x; - result->p[1].y = pt2->y; + statlseg_construct(result, pt1, pt2); PG_RETURN_LSEG_P(result); } /* like lseg_construct, but assume space already allocated */ -static void +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. @@ -2045,79 +2036,43 @@ lseg_length(PG_FUNCTION_ARGS) /* ** 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_intersect_internal(l1, l2)); + PG_RETURN_BOOL(lseg_interpt_lseg(NULL, l1, l2)); } -static bool -lseg_intersect_internal(LSEG *l1, LSEG *l2) -{ - LINE ln; - Point *interpt; - bool retval; - - line_construct_pts(&ln, &l2->p[0], &l2->p[1]); - interpt = interpt_sl(l1, &ln); - - if (interpt != NULL && on_ps_internal(interpt, l2)) - retval = true; /* interpt on l1 and l2 */ - else - retval = false; - return retval; -} Datum lseg_parallel(PG_FUNCTION_ARGS) { LSEG *l1 = PG_GETARG_LSEG_P(0); LSEG *l2 = PG_GETARG_LSEG_P(1); - PG_RETURN_BOOL(FPeq(point_sl(&l1->p[0], &l1->p[1]), - point_sl(&l2->p[0], &l2->p[1]))); + PG_RETURN_BOOL(FPeq(lseg_sl(l1), lseg_sl(l2))); } -/* lseg_perp() +/* * Determine if two line segments are perpendicular. - * - * This code did not get the correct answer for - * '((0,0),(0,1))'::lseg ?-| '((0,0),(1,0))'::lseg - * So, modified it to check explicitly for slope of vertical line - * returned by point_sl() and the results seem better. - * - thomas 1998-01-31 */ Datum lseg_perp(PG_FUNCTION_ARGS) { LSEG *l1 = PG_GETARG_LSEG_P(0); LSEG *l2 = PG_GETARG_LSEG_P(1); - double m1, - m2; - m1 = point_sl(&(l1->p[0]), &(l1->p[1])); - m2 = point_sl(&(l2->p[0]), &(l2->p[1])); - -#ifdef GEODEBUG - printf("lseg_perp- slopes are %g and %g\n", m1, m2); -#endif - if (FPzero(m1)) - PG_RETURN_BOOL(FPeq(m2, DBL_MAX)); - else if (FPzero(m2)) - PG_RETURN_BOOL(FPeq(m1, DBL_MAX)); - - PG_RETURN_BOOL(FPeq(m1 / m2, -1.0)); + 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)); } @@ -2129,36 +2084,32 @@ lseg_horizontal(PG_FUNCTION_ARGS) 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(FPeq(l1->p[0].x, l2->p[0].x) && - FPeq(l1->p[0].y, l2->p[0].y) && - FPeq(l1->p[1].x, l2->p[1].x) && - FPeq(l1->p[1].y, l2->p[1].y)); + 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(!FPeq(l1->p[0].x, l2->p[0].x) || - !FPeq(l1->p[0].y, l2->p[0].y) || - !FPeq(l1->p[1].x, l2->p[1].x) || - !FPeq(l1->p[1].y, l2->p[1].y)); + 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]))); @@ -2203,126 +2154,81 @@ lseg_ge(PG_FUNCTION_ARGS) * 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_dt(l1, l2)); -} - -/* lseg_dt() - * Distance between two line segments. - * Must check both sets of endpoints to ensure minimum distance is found. - * - thomas 1998-02-01 - */ -static double -lseg_dt(LSEG *l1, LSEG *l2) -{ - double result, - d; - - if (lseg_intersect_internal(l1, l2)) - return 0.0; - - d = dist_ps_internal(&l1->p[0], l2); - result = d; - d = dist_ps_internal(&l1->p[1], l2); - result = Min(result, d); - d = dist_ps_internal(&l2->p[0], l1); - result = Min(result, d); - d = dist_ps_internal(&l2->p[1], l1); - result = Min(result, d); - - return result; + 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 = (lseg->p[0].x + lseg->p[1].x) / 2.0; result->y = (lseg->p[0].y + lseg->p[1].y) / 2.0; PG_RETURN_POINT_P(result); } -static Point * -lseg_interpt_internal(LSEG *l1, LSEG *l2) + +/* + * Find the intersection point of two segments (if any). + * + * This returns true if two line segments intersect, false if they do not. + * This also sets the intersection point to *result, if it is not NULL. + * 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 *result; - LINE tmp1, - tmp2; + Point interpt; + LINE tmp; + + line_construct(&tmp, &l2->p[0], lseg_sl(l2)); + if (!lseg_interpt_line(&interpt, l1, &tmp)) + return false; /* - * Find the intersection of the appropriate lines, if any. + * If the line intersection point isn't within l2, there is no valid + * segment intersection point at all. */ - line_construct_pts(&tmp1, &l1->p[0], &l1->p[1]); - line_construct_pts(&tmp2, &l2->p[0], &l2->p[1]); - result = line_interpt_internal(&tmp1, &tmp2); - if (!PointerIsValid(result)) - return NULL; + if (!lseg_contain_point(l2, &interpt)) + return false; - /* - * If the line intersection point isn't within l1 (or equivalently l2), - * there is no valid segment intersection point at all. - */ - if (!on_ps_internal(result, l1) || - !on_ps_internal(result, l2)) - { - pfree(result); - return NULL; - } + if (result != NULL) + *result = interpt; - /* - * If there is an intersection, then check explicitly for matching - * endpoints since there may be rounding effects with annoying lsb - * residue. - tgl 1997-07-09 - */ - if ((FPeq(l1->p[0].x, l2->p[0].x) && FPeq(l1->p[0].y, l2->p[0].y)) || - (FPeq(l1->p[0].x, l2->p[1].x) && FPeq(l1->p[0].y, l2->p[1].y))) - { - result->x = l1->p[0].x; - result->y = l1->p[0].y; - } - else if ((FPeq(l1->p[1].x, l2->p[0].x) && FPeq(l1->p[1].y, l2->p[0].y)) || - (FPeq(l1->p[1].x, l2->p[1].x) && FPeq(l1->p[1].y, l2->p[1].y))) - { - result->x = l1->p[1].x; - result->y = l1->p[1].y; - } - - return result; + return true; } -/* lseg_interpt - - * Find the intersection point of two segments (if any). - */ Datum lseg_interpt(PG_FUNCTION_ARGS) { LSEG *l1 = PG_GETARG_LSEG_P(0); LSEG *l2 = PG_GETARG_LSEG_P(1); Point *result; - result = lseg_interpt_internal(l1, l2); - if (!PointerIsValid(result)) - PG_RETURN_NULL(); + 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. ** ***********************************************************************/ @@ -2333,215 +2239,128 @@ lseg_interpt(PG_FUNCTION_ARGS) /* * 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(dist_pl_internal(pt, line)); + PG_RETURN_FLOAT8(line_closept_point(NULL, line, pt)); } -static double -dist_pl_internal(Point *pt, LINE *line) -{ - return fabs((line->A * pt->x + line->B * pt->y + line->C) / - HYPOT(line->A, line->B)); -} /* * 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(dist_ps_internal(pt, lseg)); -} - -static double -dist_ps_internal(Point *pt, LSEG *lseg) -{ - double m; /* slope of perp. */ - LINE *ln; - double result, - tmpdist; - Point *ip; - - /* - * Construct a line perpendicular to the input segment and through the - * input point - */ - if (lseg->p[1].x == lseg->p[0].x) - m = 0; - else if (lseg->p[1].y == lseg->p[0].y) - m = (double) DBL_MAX; /* slope is infinite */ - else - m = (lseg->p[0].x - lseg->p[1].x) / (lseg->p[1].y - lseg->p[0].y); - ln = line_construct_pm(pt, m); - -#ifdef GEODEBUG - printf("dist_ps- line is A=%g B=%g C=%g from (point) slope (%f,%f) %g\n", - ln->A, ln->B, ln->C, pt->x, pt->y, m); -#endif - - /* - * Calculate distance to the line segment or to the nearest endpoint of - * the segment. - */ - - /* intersection is on the line segment? */ - if ((ip = interpt_sl(lseg, ln)) != NULL) - { - /* yes, so use distance to the intersection point */ - result = point_dt(pt, ip); -#ifdef GEODEBUG - printf("dist_ps- distance is %f to intersection point is (%f,%f)\n", - result, ip->x, ip->y); -#endif - } - else - { - /* no, so use distance to the nearer endpoint */ - result = point_dt(pt, &lseg->p[0]); - tmpdist = point_dt(pt, &lseg->p[1]); - if (tmpdist < result) - result = tmpdist; - } - - return result; + PG_RETURN_FLOAT8(lseg_closept_point(NULL, lseg, pt)); } /* * 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); float8 result = 0.0; /* keep compiler quiet */ bool have_min = false; float8 tmp; int i; LSEG lseg; - switch (path->npts) + 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++) { - case 0: - /* no points in path? then result is undefined... */ - PG_RETURN_NULL(); - case 1: - /* one point in path? then get distance between two points... */ - result = point_dt(pt, &path->p[0]); - break; - default: - /* make sure the path makes sense... */ - Assert(path->npts > 1); + int iprev; - /* - * 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 */ + } - 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 = dist_ps_internal(pt, &lseg); - if (!have_min || tmp < result) - { - result = tmp; - have_min = true; - } - } - break; + statlseg_construct(&lseg, &path->p[iprev], &path->p[i]); + tmp = lseg_closept_point(NULL, &lseg, pt); + if (!have_min || tmp < result) + { + result = tmp; + have_min = true; + } } + PG_RETURN_FLOAT8(result); } /* * 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); - float8 result; - Point *near; - near = DatumGetPointP(DirectFunctionCall2(close_pb, - PointPGetDatum(pt), - BoxPGetDatum(box))); - result = point_dt(near, pt); - - PG_RETURN_FLOAT8(result); + 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); float8 result, d2; - if (has_interpt_sl(lseg, line)) + if (lseg_interpt_line(NULL, lseg, line)) result = 0.0; else { - result = dist_pl_internal(&lseg->p[0], line); - d2 = dist_pl_internal(&lseg->p[1], line); + result = line_closept_point(NULL, line, &lseg->p[0]); + d2 = line_closept_point(NULL, line, &lseg->p[1]); /* XXX shouldn't we take the min not max? */ if (d2 > result) result = d2; } PG_RETURN_FLOAT8(result); } /* * 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); - Point *tmp; - Datum result; - tmp = DatumGetPointP(DirectFunctionCall2(close_sb, - LsegPGetDatum(lseg), - BoxPGetDatum(box))); - result = DirectFunctionCall2(dist_pb, - PointPGetDatum(tmp), - BoxPGetDatum(box)); - - PG_RETURN_DATUM(result); + 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); @@ -2577,538 +2396,509 @@ dist_cpoly(PG_FUNCTION_ARGS) } /* * 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); - float8 result; - result = dist_ppoly_internal(point, poly); - - PG_RETURN_FLOAT8(result); + 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); - float8 result; - result = dist_ppoly_internal(point, poly); - - PG_RETURN_FLOAT8(result); + PG_RETURN_FLOAT8(dist_ppoly_internal(point, poly)); } -static double +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) { #ifdef GEODEBUG printf("dist_ppoly_internal- point inside of polygon\n"); #endif 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 = dist_ps_internal(pt, &seg); + result = lseg_closept_point(NULL, &seg, pt); #ifdef GEODEBUG printf("dist_ppoly_internal- segment 0/n distance is %f\n", result); #endif /* check distances for other segments */ - for (i = 0; (i < poly->npts - 1); i++) + 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 = dist_ps_internal(pt, &seg); + d = lseg_closept_point(NULL, &seg, pt); #ifdef GEODEBUG printf("dist_ppoly_internal- segment %d distance is %f\n", (i + 1), d); #endif if (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. *-------------------------------------------------------------------*/ -/* Get intersection point of lseg and line; returns NULL if no intersection */ -static Point * -interpt_sl(LSEG *lseg, LINE *line) -{ - LINE tmp; - Point *p; - - line_construct_pts(&tmp, &lseg->p[0], &lseg->p[1]); - p = line_interpt_internal(&tmp, line); -#ifdef GEODEBUG - printf("interpt_sl- segment is (%.*g %.*g) (%.*g %.*g)\n", - DBL_DIG, lseg->p[0].x, DBL_DIG, lseg->p[0].y, DBL_DIG, lseg->p[1].x, DBL_DIG, lseg->p[1].y); - printf("interpt_sl- segment becomes line A=%.*g B=%.*g C=%.*g\n", - DBL_DIG, tmp.A, DBL_DIG, tmp.B, DBL_DIG, tmp.C); -#endif - if (PointerIsValid(p)) - { -#ifdef GEODEBUG - printf("interpt_sl- intersection point is (%.*g %.*g)\n", DBL_DIG, p->x, DBL_DIG, p->y); -#endif - if (on_ps_internal(p, lseg)) - { -#ifdef GEODEBUG - printf("interpt_sl- intersection point is on segment\n"); -#endif - } - else - p = NULL; - } - - return p; -} - -/* variant: just indicate if intersection point exists */ +/* + * Check if the line segment intersects with the line + * + * This returns true if line segment intersects with line, false if they + * do not. This also sets the intersection point to *result, if it is not + * NULL. + */ static bool -has_interpt_sl(LSEG *lseg, LINE *line) +lseg_interpt_line(Point *result, LSEG *lseg, LINE *line) { - Point *tmp; + Point interpt; + LINE tmp; - tmp = interpt_sl(lseg, line); - if (tmp) + line_construct(&tmp, &lseg->p[0], lseg_sl(lseg)); + if (!line_interpt_line(&interpt, &tmp, line)) + return false; + + if (!lseg_contain_point(lseg, &interpt)) + return false; + + if (result == NULL) return true; - return false; + + /* + * 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. *-------------------------------------------------------------------*/ -/* close_pl - +/* * The intersection point of a perpendicular of the line * through the point. + * + * This sets the closest point to the *result if it is not NULL and returns + * the distance to the closest point. */ +static float8 +line_closept_point(Point *result, LINE *line, Point *point) +{ + bool retval; + LINE tmp; + + /* Drop a perpendicular and find the intersection point */ + line_construct(&tmp, point, line_invsl(line)); + retval = line_interpt_line(result, line, &tmp); + Assert(retval); /* XXX We need something better. */ + + /* + * XXX We could use the distance to the closest point, but + * line_interpt_line() is currently giving wrong results. + */ + return fabs((line->A * point->x + line->B * point->y + line->C) / + HYPOT(line->A, line->B)); +} + Datum close_pl(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); LINE *line = PG_GETARG_LINE_P(1); Point *result; - LINE *tmp; - double invm; result = (Point *) palloc(sizeof(Point)); - if (FPzero(line->B)) /* vertical? */ - { - result->x = line->C; - result->y = pt->y; - PG_RETURN_POINT_P(result); - } - if (FPzero(line->A)) /* horizontal? */ - { - result->x = pt->x; - result->y = line->C; - PG_RETURN_POINT_P(result); - } - /* drop a perpendicular and find the intersection point */ + line_closept_point(result, line, pt); - /* invert and flip the sign on the slope to get a perpendicular */ - invm = line->B / line->A; - tmp = line_construct_pm(pt, invm); - result = line_interpt_internal(tmp, line); - Assert(result != NULL); PG_RETURN_POINT_P(result); } -/* close_ps() +/* * Closest point on line segment to specified point. - * Take the closest endpoint if the point is left, right, - * above, or below the segment, otherwise find the intersection - * point of the segment and its perpendicular through the point. * - * Some tricky code here, relying on boolean expressions - * evaluating to only zero or one to use as an array index. - * bug fixes by gthaker@atl.lmco.com; May 1, 1998 + * This sets the closest point to the *result if it is not NULL and returns + * the distance to the closest point. */ +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 = NULL; - LINE *tmp; - double invm; - int xh, - yh; + Point *result; -#ifdef GEODEBUG - printf("close_sp:pt->x %f pt->y %f\nlseg(0).x %f lseg(0).y %f lseg(1).x %f lseg(1).y %f\n", - pt->x, pt->y, lseg->p[0].x, lseg->p[0].y, - lseg->p[1].x, lseg->p[1].y); -#endif + result = (Point *) palloc(sizeof(Point)); - /* xh (or yh) is the index of upper x( or y) end point of lseg */ - /* !xh (or !yh) is the index of lower x( or y) end point of lseg */ - xh = lseg->p[0].x < lseg->p[1].x; - yh = lseg->p[0].y < lseg->p[1].y; - - if (FPeq(lseg->p[0].x, lseg->p[1].x)) /* vertical? */ - { -#ifdef GEODEBUG - printf("close_ps- segment is vertical\n"); -#endif - /* first check if point is below or above the entire lseg. */ - if (pt->y < lseg->p[!yh].y) - result = point_copy(&lseg->p[!yh]); /* below the lseg */ - else if (pt->y > lseg->p[yh].y) - result = point_copy(&lseg->p[yh]); /* above the lseg */ - if (result != NULL) - PG_RETURN_POINT_P(result); - - /* point lines along (to left or right) of the vertical lseg. */ - - result = (Point *) palloc(sizeof(Point)); - result->x = lseg->p[0].x; - result->y = pt->y; - PG_RETURN_POINT_P(result); - } - else if (FPeq(lseg->p[0].y, lseg->p[1].y)) /* horizontal? */ - { -#ifdef GEODEBUG - printf("close_ps- segment is horizontal\n"); -#endif - /* first check if point is left or right of the entire lseg. */ - if (pt->x < lseg->p[!xh].x) - result = point_copy(&lseg->p[!xh]); /* left of the lseg */ - else if (pt->x > lseg->p[xh].x) - result = point_copy(&lseg->p[xh]); /* right of the lseg */ - if (result != NULL) - PG_RETURN_POINT_P(result); - - /* point lines along (at top or below) the horiz. lseg. */ - result = (Point *) palloc(sizeof(Point)); - result->x = pt->x; - result->y = lseg->p[0].y; - PG_RETURN_POINT_P(result); - } - - /* - * vert. and horiz. cases are down, now check if the closest point is one - * of the end points or someplace on the lseg. - */ - - invm = -1.0 / point_sl(&(lseg->p[0]), &(lseg->p[1])); - tmp = line_construct_pm(&lseg->p[!yh], invm); /* lower edge of the - * "band" */ - if (pt->y < (tmp->A * pt->x + tmp->C)) - { /* we are below the lower edge */ - result = point_copy(&lseg->p[!yh]); /* below the lseg, take lower end - * pt */ -#ifdef GEODEBUG - printf("close_ps below: tmp A %f B %f C %f\n", - tmp->A, tmp->B, tmp->C); -#endif - PG_RETURN_POINT_P(result); - } - tmp = line_construct_pm(&lseg->p[yh], invm); /* upper edge of the - * "band" */ - if (pt->y > (tmp->A * pt->x + tmp->C)) - { /* we are below the lower edge */ - result = point_copy(&lseg->p[yh]); /* above the lseg, take higher end - * pt */ -#ifdef GEODEBUG - printf("close_ps above: tmp A %f B %f C %f\n", - tmp->A, tmp->B, tmp->C); -#endif - PG_RETURN_POINT_P(result); - } - - /* - * at this point the "normal" from point will hit lseg. The closest point - * will be somewhere on the lseg - */ - tmp = line_construct_pm(pt, invm); -#ifdef GEODEBUG - printf("close_ps- tmp A %f B %f C %f\n", - tmp->A, tmp->B, tmp->C); -#endif - result = interpt_sl(lseg, tmp); - - /* - * ordinarily we should always find an intersection point, but that could - * fail in the presence of NaN coordinates, and perhaps even from simple - * roundoff issues. Return a SQL NULL if so. - */ - if (result == NULL) + if (isnan(lseg_closept_point(result, lseg, pt))) PG_RETURN_NULL(); -#ifdef GEODEBUG - printf("close_ps- result.x %f result.y %f\n", result->x, result->y); -#endif PG_RETURN_POINT_P(result); } -/* close_lseg() +/* * Closest point to l1 on l2. + * + * This sets the closest point to the *result if it is not NULL and returns + * the distance to the closest point. + * + * XXX This function is wrong. If must never set the *result to a point on + * the second segment. */ +static float8 +lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2) +{ + Point point; + double dist; + double d; + + d = lseg_closept_point(NULL, l1, &l2->p[0]); + dist = d; + if (result != NULL) + *result = l2->p[0]; + + d = lseg_closept_point(NULL, l1, &l2->p[1]); + if (d < dist) + { + dist = d; + if (result != NULL) + *result = l2->p[1]; + } + + if (lseg_closept_point(&point, l2, &l1->p[0]) < dist) + d = lseg_closept_point(result, l1, &point); + + if (lseg_closept_point(&point, l2, &l1->p[1]) < dist) + d = lseg_closept_point(result, l1, &point); + + if (d < dist) + dist = d; + + return dist; +} + Datum close_lseg(PG_FUNCTION_ARGS) { LSEG *l1 = PG_GETARG_LSEG_P(0); LSEG *l2 = PG_GETARG_LSEG_P(1); - Point *result = NULL; - Point point; - double dist; - double d; + Point *result; - d = dist_ps_internal(&l1->p[0], l2); - dist = d; - memcpy(&point, &l1->p[0], sizeof(Point)); + result = (Point *) palloc(sizeof(Point)); - if ((d = dist_ps_internal(&l1->p[1], l2)) < dist) - { - dist = d; - memcpy(&point, &l1->p[1], sizeof(Point)); - } - - if (dist_ps_internal(&l2->p[0], l1) < dist) - { - result = DatumGetPointP(DirectFunctionCall2(close_ps, - PointPGetDatum(&l2->p[0]), - LsegPGetDatum(l1))); - memcpy(&point, result, sizeof(Point)); - result = DatumGetPointP(DirectFunctionCall2(close_ps, - PointPGetDatum(&point), - LsegPGetDatum(l2))); - } - - if (dist_ps_internal(&l2->p[1], l1) < dist) - { - result = DatumGetPointP(DirectFunctionCall2(close_ps, - PointPGetDatum(&l2->p[1]), - LsegPGetDatum(l1))); - memcpy(&point, result, sizeof(Point)); - result = DatumGetPointP(DirectFunctionCall2(close_ps, - PointPGetDatum(&point), - LsegPGetDatum(l2))); - } - - if (result == NULL) - result = point_copy(&point); + lseg_closept_lseg(result, l2, l1); PG_RETURN_POINT_P(result); } -/* close_pb() + +/* * Closest point on or in box to specified point. + * + * This sets the closest point to the *result if it is not NULL and returns + * the distance to the closest point. */ -Datum -close_pb(PG_FUNCTION_ARGS) +static float8 +box_closept_point(Point *result, BOX *box, Point *pt) { - Point *pt = PG_GETARG_POINT_P(0); - BOX *box = PG_GETARG_BOX_P(1); - LSEG lseg, - seg; - Point point; + LSEG lseg; + Point point, + closept; double dist, d; - if (DatumGetBool(DirectFunctionCall2(on_pb, - PointPGetDatum(pt), - BoxPGetDatum(box)))) - PG_RETURN_POINT_P(pt); + 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 = dist_ps_internal(pt, &lseg); + dist = lseg_closept_point(result, &lseg, pt); - statlseg_construct(&seg, &box->high, &point); - if ((d = dist_ps_internal(pt, &seg)) < dist) + statlseg_construct(&lseg, &box->high, &point); + d = lseg_closept_point(&closept, &lseg, pt); + if (d < dist) { dist = d; - memcpy(&lseg, &seg, sizeof(lseg)); + if (result != NULL) + *result = closept; } point.x = box->high.x; point.y = box->low.y; - statlseg_construct(&seg, &box->low, &point); - if ((d = dist_ps_internal(pt, &seg)) < dist) + statlseg_construct(&lseg, &box->low, &point); + d = lseg_closept_point(&closept, &lseg, pt); + if (d < dist) { dist = d; - memcpy(&lseg, &seg, sizeof(lseg)); + if (result != NULL) + *result = closept; } - statlseg_construct(&seg, &box->high, &point); - if ((d = dist_ps_internal(pt, &seg)) < dist) + statlseg_construct(&lseg, &box->high, &point); + d = lseg_closept_point(&closept, &lseg, pt); + if (d < dist) { dist = d; - memcpy(&lseg, &seg, sizeof(lseg)); + if (result != NULL) + *result = closept; } - PG_RETURN_DATUM(DirectFunctionCall2(close_ps, - PointPGetDatum(pt), - LsegPGetDatum(&lseg))); + 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)); + + box_closept_point(result, box, pt); + + 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 = interpt_sl(lseg, line); - if (result) + result = (Point *) palloc(sizeof(Point)); + + if (lseg_interpt_line(result, lseg, line)) PG_RETURN_POINT_P(result); - d1 = dist_pl_internal(&lseg->p[0], line); - d2 = dist_pl_internal(&lseg->p[1], line); + d1 = line_closept_point(NULL, line, &lseg->p[0]); + d2 = line_closept_point(NULL, line, &lseg->p[1]); if (d1 < d2) - result = point_copy(&lseg->p[0]); + *result = lseg->p[0]; else - result = point_copy(&lseg->p[1]); + *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(); } -/* close_ls() +/* * Closest point on line segment to line. + * + * This sets the closest point to the *result if it is not NULL and returns + * the distance to the closest point. + * + * NOTE When the lines are parallel, endpoints of one of the line segment + * are FPeq(), in presence of NaN or Infinitive 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; - float8 d1, - d2; - result = interpt_sl(lseg, line); - if (result) - PG_RETURN_POINT_P(result); + result = (Point *) palloc(sizeof(Point)); - d1 = dist_pl_internal(&lseg->p[0], line); - d2 = dist_pl_internal(&lseg->p[1], line); - if (d1 < d2) - result = point_copy(&lseg->p[0]); - else - result = point_copy(&lseg->p[1]); + lseg_closept_line(result, lseg, line); PG_RETURN_POINT_P(result); } -/* close_sb() + +/* * Closest point on or in box to line segment. + * + * This sets the closest point to the *result if it is not NULL and returns + * the distance to the closest point. */ -Datum -close_sb(PG_FUNCTION_ARGS) +static float8 +box_closept_lseg(Point *result, BOX *box, LSEG *lseg) { - LSEG *lseg = PG_GETARG_LSEG_P(0); - BOX *box = PG_GETARG_BOX_P(1); - Point point; - LSEG bseg, - seg; + Point point, + closept; + LSEG bseg; double dist, d; - /* segment intersects box? then just return closest point to center */ - if (DatumGetBool(DirectFunctionCall2(inter_sb, - LsegPGetDatum(lseg), - BoxPGetDatum(box)))) - { - box_cn(&point, box); - PG_RETURN_DATUM(DirectFunctionCall2(close_ps, - PointPGetDatum(&point), - LsegPGetDatum(lseg))); - } + 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_dt(lseg, &bseg); + dist = lseg_closept_lseg(result, &bseg, lseg); - statlseg_construct(&seg, &box->high, &point); - if ((d = lseg_dt(lseg, &seg)) < dist) + statlseg_construct(&bseg, &box->high, &point); + d = lseg_closept_lseg(&closept, &bseg, lseg); + if (d < dist) { dist = d; - memcpy(&bseg, &seg, sizeof(bseg)); + if (result != NULL) + *result = closept; } point.x = box->high.x; point.y = box->low.y; - statlseg_construct(&seg, &box->low, &point); - if ((d = lseg_dt(lseg, &seg)) < dist) + statlseg_construct(&bseg, &box->low, &point); + d = lseg_closept_lseg(&closept, &bseg, lseg); + if (d < dist) { dist = d; - memcpy(&bseg, &seg, sizeof(bseg)); + if (result != NULL) + *result = closept; } - statlseg_construct(&seg, &box->high, &point); - if ((d = lseg_dt(lseg, &seg)) < dist) + statlseg_construct(&bseg, &box->high, &point); + d = lseg_closept_lseg(&closept, &bseg, lseg); + if (d < dist) { dist = d; - memcpy(&bseg, &seg, sizeof(bseg)); + if (result != NULL) + *result = closept; } - /* OK, we now have the closest line segment on the box boundary */ - PG_RETURN_DATUM(DirectFunctionCall2(close_lseg, - LsegPGetDatum(lseg), - LsegPGetDatum(&bseg))); + 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)); + + box_closept_lseg(result, box, lseg); + + 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, @@ -3116,71 +2906,87 @@ close_lb(PG_FUNCTION_ARGS) errmsg("function \"close_lb\" not implemented"))); PG_RETURN_NULL(); } /*--------------------------------------------------------------------- * on_ * Whether one object lies completely within another. *-------------------------------------------------------------------*/ -/* on_pl - +/* * Does the point satisfy the equation? */ +static bool +line_contain_point(LINE *line, Point *point) +{ + return FPzero(line->A * point->x + 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(FPzero(line->A * pt->x + line->B * pt->y + line->C)); + PG_RETURN_BOOL(line_contain_point(line, pt)); } -/* on_ps - +/* * 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(on_ps_internal(pt, lseg)); + PG_RETURN_BOOL(lseg_contain_point(lseg, pt)); } + +/* + * Check whether the point is in the box or on its border + */ static bool -on_ps_internal(Point *pt, LSEG *lseg) +box_contain_point(BOX *box, Point *point) { - return FPeq(point_dt(pt, &lseg->p[0]) + point_dt(pt, &lseg->p[1]), - point_dt(&lseg->p[0], &lseg->p[1])); + 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(pt->x <= box->high.x && pt->x >= box->low.x && - pt->y <= box->high.y && pt->y >= box->low.y); + 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(pt->x <= box->high.x && pt->x >= box->low.x && - pt->y <= box->high.y && pt->y >= box->low.y); + 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 @@ -3210,158 +3016,184 @@ on_ppath(PG_FUNCTION_ARGS) 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(DatumGetBool(DirectFunctionCall2(on_pl, - PointPGetDatum(&lseg->p[0]), - LinePGetDatum(line))) && - DatumGetBool(DirectFunctionCall2(on_pl, - PointPGetDatum(&lseg->p[1]), - LinePGetDatum(line)))); + 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(DatumGetBool(DirectFunctionCall2(on_pb, - PointPGetDatum(&lseg->p[0]), - BoxPGetDatum(box))) && - DatumGetBool(DirectFunctionCall2(on_pb, - PointPGetDatum(&lseg->p[1]), - BoxPGetDatum(box)))); + 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(has_interpt_sl(lseg, line)); + PG_RETURN_BOOL(lseg_interpt_line(NULL, lseg, line)); } -/* inter_sb() + +/* * 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 */ -Datum -inter_sb(PG_FUNCTION_ARGS) +static bool +box_interpt_lseg(Point *result, BOX *box, LSEG *lseg) { - LSEG *lseg = PG_GETARG_LSEG_P(0); - BOX *box = PG_GETARG_BOX_P(1); BOX lbox; LSEG bseg; Point point; lbox.low.x = Min(lseg->p[0].x, lseg->p[1].x); lbox.low.y = Min(lseg->p[0].y, lseg->p[1].y); lbox.high.x = Max(lseg->p[0].x, lseg->p[1].x); lbox.high.y = Max(lseg->p[0].y, lseg->p[1].y); /* nothing close to overlap? then not going to intersect */ if (!box_ov(&lbox, box)) - PG_RETURN_BOOL(false); + 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 (DatumGetBool(DirectFunctionCall2(on_pb, - PointPGetDatum(&lseg->p[0]), - BoxPGetDatum(box))) || - DatumGetBool(DirectFunctionCall2(on_pb, - PointPGetDatum(&lseg->p[1]), - BoxPGetDatum(box)))) - PG_RETURN_BOOL(true); + 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_intersect_internal(&bseg, lseg)) - PG_RETURN_BOOL(true); + if (lseg_interpt_lseg(NULL, &bseg, lseg)) + return true; statlseg_construct(&bseg, &box->high, &point); - if (lseg_intersect_internal(&bseg, lseg)) - PG_RETURN_BOOL(true); + 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_intersect_internal(&bseg, lseg)) - PG_RETURN_BOOL(true); + if (lseg_interpt_lseg(NULL, &bseg, lseg)) + return true; statlseg_construct(&bseg, &box->high, &point); - if (lseg_intersect_internal(&bseg, lseg)) - PG_RETURN_BOOL(true); + if (lseg_interpt_lseg(NULL, &bseg, lseg)) + return true; /* if we dropped through, no two segs intersected */ - PG_RETURN_BOOL(false); + 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 (has_interpt_sl(&bseg, line)) + 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 (has_interpt_sl(&bseg, line)) + 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 (has_interpt_sl(&bseg, line)) + 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 (has_interpt_sl(&bseg, line)) + 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. @@ -3374,42 +3206,40 @@ inter_lb(PG_FUNCTION_ARGS) *---------------------------------------------------------------------*/ static void make_bound_box(POLYGON *poly) { int i; double x1, y1, x2, y2; - if (poly->npts > 0) - { - x2 = x1 = poly->p[0].x; - y2 = y1 = poly->p[0].y; - for (i = 1; i < poly->npts; i++) - { - if (poly->p[i].x < x1) - x1 = poly->p[i].x; - if (poly->p[i].x > x2) - x2 = poly->p[i].x; - if (poly->p[i].y < y1) - y1 = poly->p[i].y; - if (poly->p[i].y > y2) - y2 = poly->p[i].y; - } + Assert(poly->npts > 0); - box_fill(&(poly->boundbox), x1, x2, y1, y2); + x1 = x2 = poly->p[0].x; + y2 = y1 = poly->p[0].y; + for (i = 1; i < poly->npts; i++) + { + if (poly->p[i].x < x1) + x1 = poly->p[i].x; + if (poly->p[i].x > x2) + x2 = poly->p[i].x; + if (poly->p[i].y < y1) + y1 = poly->p[i].y; + if (poly->p[i].y > y2) + y2 = poly->p[i].y; } - else - ereport(ERROR, - (errcode(ERRCODE_INVALID_PARAMETER_VALUE), - errmsg("cannot create bounding box for empty polygon"))); + + 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)" *------------------------------------------------------------------*/ @@ -3739,65 +3569,65 @@ poly_same(PG_FUNCTION_ARGS) /*----------------------------------------------------------------- * 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 = (polya->npts > 0 && polyb->npts > 0 && - box_ov(&polya->boundbox, &polyb->boundbox)) ? true : false; + 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 == false; ia++) + 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 == false; ib++) + for (ib = 0; ib < polyb->npts && !result; ib++) { sb.p[1] = polyb->p[ib]; - result = lseg_intersect_internal(&sa, &sb); + 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 == false) + if (!result) { - result = (point_inside(polya->p, polyb->npts, polyb->p) - || + 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); @@ -3817,36 +3647,35 @@ poly_overlap(PG_FUNCTION_ARGS) 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; -#define POINTEQ(pt1, pt2) (FPeq((pt1)->x, (pt2)->x) && FPeq((pt1)->y, (pt2)->y)) - if (POINTEQ(a, s->p)) + if (point_eq_point(a, s->p)) { - if (on_ps_internal(s->p + 1, &t)) + if (lseg_contain_point(&t, s->p + 1)) return lseg_inside_poly(b, s->p + 1, poly, start); } - else if (POINTEQ(a, s->p + 1)) + else if (point_eq_point(a, s->p + 1)) { - if (on_ps_internal(s->p, &t)) + if (lseg_contain_point(&t, s->p)) return lseg_inside_poly(b, s->p, poly, start); } - else if (on_ps_internal(s->p, &t)) + else if (lseg_contain_point(&t, s->p)) { return lseg_inside_poly(b, s->p, poly, start); } - else if (on_ps_internal(s->p + 1, &t)) + 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 @@ -3860,50 +3689,49 @@ lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start) 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; + Point interpt; CHECK_FOR_INTERRUPTS(); s.p[1] = poly->p[i]; - if (on_ps_internal(t.p, &s)) + if (lseg_contain_point(&s, t.p)) { - if (on_ps_internal(t.p + 1, &s)) + 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 (on_ps_internal(t.p + 1, &s)) + 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 ((interpt = lseg_interpt_internal(&t, &s)) != NULL) + 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); + res = lseg_inside_poly(t.p, &interpt, poly, i + 1); if (res) - res = lseg_inside_poly(t.p + 1, interpt, poly, i + 1); - pfree(interpt); + res = lseg_inside_poly(t.p + 1, &interpt, poly, i + 1); } s.p[0] = s.p[1]; } if (res && !intersection) { Point p; /* @@ -3915,74 +3743,86 @@ lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start) res = point_inside(&p, poly->npts, poly->p); } return res; } /*----------------------------------------------------------------- * Determine if polygon A contains polygon B. *-----------------------------------------------------------------*/ +static bool +poly_contain_poly(POLYGON *polya, POLYGON *polyb) +{ + int i; + LSEG s; + + Assert(polya->npts > 0 && polyb->npts > 0); + + /* + * Quick check to see if bounding box is contained. + */ + if (!box_contain_box(&polya->boundbox, &polyb->boundbox)) + return false; + + s.p[0] = polyb->p[polyb->npts - 1]; + + for (i = 0; i < polyb->npts; i++) + { + s.p[1] = polyb->p[i]; + if (!lseg_inside_poly(s.p, s.p + 1, polya, 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; - /* - * Quick check to see if bounding box is contained. - */ - if (polya->npts > 0 && polyb->npts > 0 && - DatumGetBool(DirectFunctionCall2(box_contain, - BoxPGetDatum(&polya->boundbox), - BoxPGetDatum(&polyb->boundbox)))) - { - int i; - LSEG s; - - s.p[0] = polyb->p[polyb->npts - 1]; - result = true; - - for (i = 0; i < polyb->npts && result; i++) - { - s.p[1] = polyb->p[i]; - result = lseg_inside_poly(s.p, s.p + 1, polya, 0); - s.p[0] = s.p[1]; - } - } - else - { - result = false; - } + 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) { - Datum polya = PG_GETARG_DATUM(0); - Datum polyb = PG_GETARG_DATUM(1); + 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 */ - PG_RETURN_DATUM(DirectFunctionCall2(poly_contain, polyb, polya)); + 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); @@ -4018,170 +3858,215 @@ poly_distance(PG_FUNCTION_ARGS) ** ** Routines for 2D points. ** ***********************************************************************/ Datum construct_point(PG_FUNCTION_ARGS) { float8 x = PG_GETARG_FLOAT8(0); float8 y = PG_GETARG_FLOAT8(1); + Point *result; - PG_RETURN_POINT_P(point_construct(x, y)); + 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, + pt1->x + pt2->x, + 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)); - result->x = (p1->x + p2->x); - result->y = (p1->y + p2->y); + 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, + pt1->x - pt2->x, + 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)); - result->x = (p1->x - p2->x); - result->y = (p1->y - p2->y); + 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, + (pt1->x * pt2->x) - (pt1->y * pt2->y), + (pt1->x * pt2->y) + (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)); - result->x = (p1->x * p2->x) - (p1->y * p2->y); - result->y = (p1->x * p2->y) + (p1->y * p2->x); + point_mul_point(result, p1, p2); PG_RETURN_POINT_P(result); } + +static inline void +point_div_point(Point *result, Point *pt1, Point *pt2) +{ + double div; + + div = (pt2->x * pt2->x) + (pt2->y * pt2->y); + + if (div == 0.0) + ereport(ERROR, + (errcode(ERRCODE_DIVISION_BY_ZERO), + errmsg("division by zero"))); + + point_construct(result, + ((pt1->x * pt2->x) + (pt1->y * pt2->y)) / div, + ((pt2->x * pt1->y) - (pt2->y * pt1->x)) / div); +} + Datum point_div(PG_FUNCTION_ARGS) { Point *p1 = PG_GETARG_POINT_P(0); Point *p2 = PG_GETARG_POINT_P(1); Point *result; - double div; result = (Point *) palloc(sizeof(Point)); - div = (p2->x * p2->x) + (p2->y * p2->y); - - if (div == 0.0) - ereport(ERROR, - (errcode(ERRCODE_DIVISION_BY_ZERO), - errmsg("division by zero"))); - - result->x = ((p1->x * p2->x) + (p1->y * p2->y)) / div; - result->y = ((p2->x * p1->y) - (p2->y * p1->x)) / div; + 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; - PG_RETURN_BOX_P(box_construct(p1->x, p2->x, p1->y, p2->y)); + 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; - PG_RETURN_BOX_P(box_construct((box->high.x + p->x), - (box->low.x + p->x), - (box->high.y + p->y), - (box->low.y + p->y))); + 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; - PG_RETURN_BOX_P(box_construct((box->high.x - p->x), - (box->low.x - p->x), - (box->high.y - p->y), - (box->low.y - p->y))); + 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; + Point high, + low; - high = DatumGetPointP(DirectFunctionCall2(point_mul, - PointPGetDatum(&box->high), - PointPGetDatum(p))); - low = DatumGetPointP(DirectFunctionCall2(point_mul, - PointPGetDatum(&box->low), - PointPGetDatum(p))); + result = (BOX *) palloc(sizeof(BOX)); - result = box_construct(high->x, low->x, high->y, low->y); + 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; + Point high, + low; - high = DatumGetPointP(DirectFunctionCall2(point_div, - PointPGetDatum(&box->high), - PointPGetDatum(p))); - low = DatumGetPointP(DirectFunctionCall2(point_div, - PointPGetDatum(&box->low), - PointPGetDatum(p))); + result = (BOX *) palloc(sizeof(BOX)); - result = box_construct(high->x, low->x, high->y, low->y); + 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) { @@ -4277,83 +4162,63 @@ path_add(PG_FUNCTION_ARGS) * 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++) - { - path->p[i].x += point->x; - path->p[i].y += point->y; - } + 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++) - { - path->p[i].x -= point->x; - path->p[i].y -= point->y; - } + 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); - Point *p; int i; for (i = 0; i < path->npts; i++) - { - p = DatumGetPointP(DirectFunctionCall2(point_mul, - PointPGetDatum(&path->p[i]), - PointPGetDatum(point))); - path->p[i].x = p->x; - path->p[i].y = p->y; - } + 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); - Point *p; int i; for (i = 0; i < path->npts; i++) - { - p = DatumGetPointP(DirectFunctionCall2(point_div, - PointPGetDatum(&path->p[i]), - PointPGetDatum(point))); - path->p[i].x = p->x; - path->p[i].y = p->y; - } + 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); @@ -4414,42 +4279,40 @@ 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); - Datum result; - CIRCLE *circle; + Point *result; + CIRCLE circle; - circle = DatumGetCircleP(DirectFunctionCall1(poly_circle, - PolygonPGetDatum(poly))); - result = DirectFunctionCall1(circle_center, - CirclePGetDatum(circle)); + result = (Point *) palloc(sizeof(Point)); - PG_RETURN_DATUM(result); + 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; - if (poly->npts < 1) - PG_RETURN_NULL(); - - box = box_copy(&poly->boundbox); + 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) @@ -4467,22 +4330,21 @@ box_poly(PG_FUNCTION_ARGS) 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_fill(&poly->boundbox, box->high.x, box->low.x, - box->high.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; @@ -4557,22 +4419,21 @@ circle_in(PG_FUNCTION_ARGS) circle->radius = single_decode(s, &s, "circle", str); if (circle->radius < 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))) + 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\"", @@ -4656,22 +4517,21 @@ circle_send(PG_FUNCTION_ARGS) /* circles identical? */ Datum circle_same(PG_FUNCTION_ARGS) { CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); PG_RETURN_BOOL(FPeq(circle1->radius, circle2->radius) && - FPeq(circle1->center.x, circle2->center.x) && - FPeq(circle1->center.y, circle2->center.y)); + 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); @@ -4730,32 +4590,34 @@ circle_overright(PG_FUNCTION_ARGS) } /* 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) + circle1->radius), circle2->radius)); + PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center), + 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) + circle2->radius), circle1->radius)); + PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center), + 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); @@ -4858,107 +4720,83 @@ circle_ge(PG_FUNCTION_ARGS) CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); PG_RETURN_BOOL(FPge(circle_ar(circle1), circle_ar(circle2))); } /*---------------------------------------------------------- * "Arithmetic" operators on circles. *---------------------------------------------------------*/ -static CIRCLE * -circle_copy(CIRCLE *circle) -{ - CIRCLE *result; - - if (!PointerIsValid(circle)) - return NULL; - - result = (CIRCLE *) palloc(sizeof(CIRCLE)); - memcpy((char *) result, (char *) circle, sizeof(CIRCLE)); - return result; -} - - /* 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_copy(circle); + result = (CIRCLE *) palloc(sizeof(CIRCLE)); - result->center.x += point->x; - result->center.y += point->y; + 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_copy(circle); + result = (CIRCLE *) palloc(sizeof(CIRCLE)); - result->center.x -= point->x; - result->center.y -= point->y; + 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; - Point *p; - result = circle_copy(circle); + result = (CIRCLE *) palloc(sizeof(CIRCLE)); - p = DatumGetPointP(DirectFunctionCall2(point_mul, - PointPGetDatum(&circle->center), - PointPGetDatum(point))); - result->center.x = p->x; - result->center.y = p->y; - result->radius *= HYPOT(point->x, point->y); + point_mul_point(&result->center, &circle->center, point); + result->radius = 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; - Point *p; - result = circle_copy(circle); + result = (CIRCLE *) palloc(sizeof(CIRCLE)); - p = DatumGetPointP(DirectFunctionCall2(point_div, - PointPGetDatum(&circle->center), - PointPGetDatum(point))); - result->center.x = p->x; - result->center.y = p->y; - result->radius /= HYPOT(point->x, point->y); + point_div_point(&result->center, &circle->center, point); + result->radius = 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) { @@ -4993,22 +4831,22 @@ circle_radius(PG_FUNCTION_ARGS) /* 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 = point_dt(&circle1->center, &circle2->center) - - (circle1->radius + circle2->radius); + result = point_dt(&circle1->center, &circle2->center) - + (circle1->radius + circle2->radius); if (result < 0) result = 0; PG_RETURN_FLOAT8(result); } Datum circle_contain_pt(PG_FUNCTION_ARGS) { CIRCLE *circle = PG_GETARG_CIRCLE_P(0); @@ -5190,56 +5028,60 @@ circle_poly(PG_FUNCTION_ARGS) angle = i * anglestep; poly->p[i].x = circle->center.x - (circle->radius * cos(angle)); poly->p[i].y = circle->center.y + (circle->radius * sin(angle)); } make_bound_box(poly); PG_RETURN_POLYGON_P(poly); } -/* poly_circle - convert polygon to circle +/* + * 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 /= poly->npts; + result->center.y /= poly->npts; + + for (i = 0; i < poly->npts; i++) + result->radius += point_dt(&poly->p[i], &result->center); + result->radius /= poly->npts; +} + Datum poly_circle(PG_FUNCTION_ARGS) { POLYGON *poly = PG_GETARG_POLYGON_P(0); - CIRCLE *circle; - int i; + CIRCLE *result; - if (poly->npts < 1) - ereport(ERROR, - (errcode(ERRCODE_INVALID_PARAMETER_VALUE), - errmsg("cannot convert empty polygon to circle"))); + result = (CIRCLE *) palloc(sizeof(CIRCLE)); - circle = (CIRCLE *) palloc(sizeof(CIRCLE)); + poly_to_circle(result, poly); - circle->center.x = 0; - circle->center.y = 0; - circle->radius = 0; - - for (i = 0; i < poly->npts; i++) - { - circle->center.x += poly->p[i].x; - circle->center.y += poly->p[i].y; - } - circle->center.x /= poly->npts; - circle->center.y /= poly->npts; - - for (i = 0; i < poly->npts; i++) - circle->radius += point_dt(&poly->p[i], &circle->center); - circle->radius /= poly->npts; - - PG_RETURN_CIRCLE_P(circle); + PG_RETURN_CIRCLE_P(result); } /*********************************************************************** ** ** Private routines for multiple types. ** ***********************************************************************/ /* @@ -5261,22 +5103,21 @@ point_inside(Point *p, int npts, Point *plist) double x0, y0; double prev_x, prev_y; int i = 0; double x, y; int cross, total_cross = 0; - if (npts <= 0) - return 0; + Assert(npts > 0); /* compute first polygon point relative to single point */ x0 = plist[0].x - p->x; y0 = 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++) { @@ -5371,51 +5212,48 @@ lseg_crossing(double x, double y, double prev_x, double prev_y) 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 ((FPeq(p2[i].x, p1[0].x)) - && (FPeq(p2[i].y, p1[0].y))) + 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 ((!FPeq(p2[j].x, p1[ii].x)) - || (!FPeq(p2[j].y, p1[ii].y))) + if (!point_eq_point(&p2[j], &p1[ii])) { #ifdef GEODEBUG printf("plist_same- %d failed forward match with %d\n", j, ii); #endif break; } } #ifdef GEODEBUG printf("plist_same- ii = %d/%d after forward match\n", ii, npts); #endif 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 ((!FPeq(p2[j].x, p1[ii].x)) - || (!FPeq(p2[j].y, p1[ii].y))) + if (!point_eq_point(&p2[j], &p1[ii])) { #ifdef GEODEBUG printf("plist_same- %d failed reverse match with %d\n", j, ii); #endif break; } } #ifdef GEODEBUG printf("plist_same- ii = %d/%d after reverse match\n", ii, npts); #endif diff --git a/src/backend/utils/adt/geo_spgist.c b/src/backend/utils/adt/geo_spgist.c index 3f1a755cbb..0fe40499e0 100644 --- a/src/backend/utils/adt/geo_spgist.c +++ b/src/backend/utils/adt/geo_spgist.c @@ -786,14 +786,15 @@ spg_bbox_quad_config(PG_FUNCTION_ARGS) /* * SP-GiST compress function for polygons */ Datum spg_poly_quad_compress(PG_FUNCTION_ARGS) { POLYGON *polygon = PG_GETARG_POLYGON_P(0); BOX *box; - box = box_copy(&polygon->boundbox); + box = (BOX *) palloc(sizeof(BOX)); + *box = polygon->boundbox; PG_RETURN_BOX_P(box); } diff --git a/src/include/utils/geo_decls.h b/src/include/utils/geo_decls.h index a589e4239f..0e066894cd 100644 --- a/src/include/utils/geo_decls.h +++ b/src/include/utils/geo_decls.h @@ -171,17 +171,13 @@ typedef struct #define DatumGetCircleP(X) ((CIRCLE *) DatumGetPointer(X)) #define CirclePGetDatum(X) PointerGetDatum(X) #define PG_GETARG_CIRCLE_P(n) DatumGetCircleP(PG_GETARG_DATUM(n)) #define PG_RETURN_CIRCLE_P(x) return CirclePGetDatum(x) /* * in geo_ops.c */ -/* private routines */ -extern double point_dt(Point *pt1, Point *pt2); -extern double point_sl(Point *pt1, Point *pt2); extern double pg_hypot(double x, double y); -extern BOX *box_copy(BOX *box); #endif /* GEO_DECLS_H */ diff --git a/src/test/regress/regress.c b/src/test/regress/regress.c index e14322c798..990a84eedc 100644 --- a/src/test/regress/regress.c +++ b/src/test/regress/regress.c @@ -171,22 +171,27 @@ widget_out(PG_FUNCTION_ARGS) PG_RETURN_CSTRING(str); } PG_FUNCTION_INFO_V1(pt_in_widget); Datum pt_in_widget(PG_FUNCTION_ARGS) { Point *point = PG_GETARG_POINT_P(0); WIDGET *widget = (WIDGET *) PG_GETARG_POINTER(1); + float8 distance; - PG_RETURN_BOOL(point_dt(point, &widget->center) < widget->radius); + distance = DatumGetFloat8(DirectFunctionCall2(point_distance, + PointPGetDatum(point), + PointPGetDatum(&widget->center))); + + PG_RETURN_BOOL(distance < widget->radius); } PG_FUNCTION_INFO_V1(reverse_name); Datum reverse_name(PG_FUNCTION_ARGS) { char *string = PG_GETARG_CSTRING(0); int i; int len; -- 2.14.3 (Apple Git-98)