From 745e4f59d00cd6fd7ab304c72843dc81cae498b1 Mon Sep 17 00:00:00 2001 From: Emre Hasegeli Date: Sat, 27 May 2017 11:27:18 -0400 Subject: [PATCH 1/4] geo-funcs-v04 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. We can do better than that, but it would require touching many more lines. The changes can be summarised as: * Re-use more functions to implement others * Unify *_construct and *_interpt_internal functions to obtain pre-allocated memory * Remove private functions from geo_decls.h * Switch using C99 hypot() as the comment suggested * Add comments to describe some functions --- src/backend/utils/adt/geo_ops.c | 1036 +++++++++++++++++---------------------- src/include/utils/geo_decls.h | 13 +- src/test/regress/regress.c | 11 +- 3 files changed, 466 insertions(+), 594 deletions(-) diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c index 0348855b11..fe4155d0b2 100644 --- a/src/backend/utils/adt/geo_ops.c +++ b/src/backend/utils/adt/geo_ops.c @@ -31,60 +31,77 @@ /* * Internal routines */ enum path_delim { PATH_NONE, PATH_OPEN, PATH_CLOSED }; +/* Routines for float8 */ +static inline float8 float8_slope(float8 x1, float8 x2, float8 y1, float8 y2); + +/* Routines for two-dimensional points */ +static inline void point_construct(Point *result, float8 x, float8 y); +static inline void point_add_internal(Point *result, Point *pt1, Point *pt2); +static inline void point_sub_internal(Point *result, Point *pt1, Point *pt2); +static inline void point_mul_internal(Point *result, Point *pt1, Point *pt2); +static inline void point_div_internal(Point *result, Point *pt1, Point *pt2); +static inline bool point_eq_internal(Point *pt1, Point *pt2); +static float8 point_dt(Point *pt1, Point *pt2); static int point_inside(Point *p, int npts, Point *plist); -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_copy(BOX *box); -static BOX *box_fill(BOX *result, double x1, double x2, double y1, double y2); + +/* Routines for two-dimensional 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); + +/* Routines for two-dimensional 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); + +/* Routines for two-dimensional lines */ +static inline void line_construct_pm(LINE *result, Point *pt, float8 m); +static inline void line_construct_pts(LINE *result, Point *pt1, Point *pt2); +static bool line_interpt_internal(Point *result, LINE *l1, LINE *l2); +static inline float8 line_calculate_point(LINE *line, Point *pt); + +/* Routines for two-dimensional line segments */ +static inline void statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2); +static bool lseg_interpt_line_internal(Point *result, LSEG *lseg, LINE *line); +static bool lseg_interpt_internal(Point *result, LSEG *l1, LSEG *l2); +static bool lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start); +static float8 lseg_dt(LSEG *l1, LSEG *l2); +static int lseg_crossing(double x, double y, double px, double py); static bool on_ps_internal(Point *pt, LSEG *lseg); + +/* Routines for two-dimensional polygons */ static void make_bound_box(POLYGON *poly); + +/* Unsorted */ 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 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 '(' @@ -327,20 +344,44 @@ pair_count(char *s, char delim) { ndelim++; s++; } return (ndelim % 2) ? ((ndelim + 1) / 2) : -1; } /*********************************************************************** ** + ** Routines for float8 + ** + ***********************************************************************/ + +/* + * Return slope of two points + * + * This function accepts x and y coordinates separately to let it be used + * to calculate inverse slope. To achieve that pass the values in + * (y1, y2, x2, x1) order. + */ +static inline float8 +float8_slope(float8 x1, float8 x2, float8 y1, float8 y2) +{ + if (FPeq(x1, x2)) + return DBL_MAX; + if (FPeq(y1, y2)) + return 0.0; + return (y1 - y2) / (x1 - x2); +} + + +/*********************************************************************** + ** ** Routines for two-dimensional boxes. ** ***********************************************************************/ /*---------------------------------------------------------- * Formatting and conversion routines. *---------------------------------------------------------*/ /* box_in - convert a string to internal form. * @@ -434,89 +475,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 - */ -static 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_internal(&box1->high, &box2->high) && + point_eq_internal(&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); @@ -751,51 +764,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)); @@ -934,22 +947,22 @@ line_in(PG_FUNCTION_ARGS) (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_internal(&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]); } PG_RETURN_LINE_P(line); } @@ -997,128 +1010,104 @@ 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_pm(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) +static inline void +line_construct_pts(LINE *result, 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 - } + float8 m; + + m = float8_slope(pt1->x, pt2->x, pt1->y, pt2->y); + line_construct_pm(result, pt1, m); } -/* line_construct_pp() - * two points + +/* + * Construct a line from 2 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); PG_RETURN_LINE_P(result); } +/* + * Calculate the line equation for a point + * + * This returns the result of the line equation Ax + By + C. The result + * needs to be 0 for the point to be on the line. + */ +static inline float8 +line_calculate_point(LINE *line, Point *pt) +{ + return line->A * pt->x + line->B * pt->y + line->C; +} + + /*---------------------------------------------------------- * 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_internal(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_internal(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)); @@ -1172,96 +1161,94 @@ line_eq(PG_FUNCTION_ARGS) /* 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_internal(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 = dist_pl_internal(&tmp, l2); 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_internal(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_internal(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 { - x = (l1->C - l2->C) / (l2->A - l1->A); + if (FPeq(l2->A, l1->A * (l2->B / l1->B))) + return false; + + if (FPzero(l2->B)) + x = l2->C; + else + 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 @@ -1609,21 +1596,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_internal(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 @@ -1664,23 +1651,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_dt(&seg1, &seg2); if (!have_min || tmp < min) { min = tmp; have_min = true; } } } if (!have_min) PG_RETURN_NULL(); @@ -1774,44 +1759,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). @@ -1870,71 +1842,74 @@ 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_internal(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_internal(pt1, pt2)); } +static inline bool +point_eq_internal(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 float8 point_dt(Point *pt1, Point *pt2) { + float8 result; + + result = hypot(pt1->x - pt2->x, pt1->y - pt2->y); + #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)); + pt1->x, pt1->y, pt2->x, pt2->y, result); #endif - return HYPOT(pt1->x - pt2->x, pt1->y - pt2->y); + + return result; } 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 -point_sl(Point *pt1, Point *pt2) -{ - return (FPeq(pt1->x, pt2->x) - ? (double) DBL_MAX - : (pt1->y - pt2->y) / (pt1->x - pt2->x)); + PG_RETURN_FLOAT8(float8_slope(pt1->x, pt2->x, pt1->y, pt2->y)); } /*********************************************************************** ** ** Routines for 2D line segments. ** ***********************************************************************/ /*---------------------------------------------------------- @@ -1946,31 +1921,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; @@ -2015,21 +1990,21 @@ lseg_construct(PG_FUNCTION_ARGS) result->p[0].x = pt1->x; result->p[0].y = pt1->y; result->p[1].x = pt2->x; result->p[1].y = pt2->y; 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; } Datum lseg_length(PG_FUNCTION_ARGS) @@ -2046,69 +2021,57 @@ 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_internal(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); + float8 m1, + m2; - PG_RETURN_BOOL(FPeq(point_sl(&l1->p[0], &l1->p[1]), - point_sl(&l2->p[0], &l2->p[1]))); + m1 = float8_slope(l1->p[0].x, l1->p[1].x, l1->p[0].y, l1->p[1].y); + m2 = float8_slope(l2->p[0].x, l2->p[1].x, l2->p[0].y, l2->p[1].y); + + PG_RETURN_BOOL(FPeq(m1, m2)); } /* 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. + * returned by float8_slope() 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])); + m1 = float8_slope(l1->p[0].x, l1->p[1].x, l1->p[0].y, l1->p[1].y); + m2 = float8_slope(l2->p[0].x, l2->p[1].x, l2->p[0].y, l2->p[1].y); #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)); @@ -2130,36 +2093,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_internal(&l1->p[0], &l2->p[0]) && + point_eq_internal(&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_internal(&l1->p[0], &l2->p[0]) || + !point_eq_internal(&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]))); @@ -2218,21 +2177,21 @@ lseg_distance(PG_FUNCTION_ARGS) * 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)) + if (lseg_interpt_internal(NULL, 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); @@ -2248,82 +2207,86 @@ lseg_center(PG_FUNCTION_ARGS) 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) +static bool +lseg_interpt_internal(Point *result, LSEG *l1, LSEG *l2) { - Point *result; + Point interpt; LINE tmp1, tmp2; /* * Find the intersection of the appropriate lines, if any. */ 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 (!line_interpt_internal(&interpt, &tmp1, &tmp2)) + 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 (!on_ps_internal(&interpt, l1) || + !on_ps_internal(&interpt, l2)) + return false; + + if (result == NULL) + return true; /* * 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; } + else + { + result->x = interpt.x; + result->y = interpt.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_internal(result, l1, l2)) + PG_RETURN_NULL(); PG_RETURN_POINT_P(result); } /*********************************************************************** ** ** Routines for position comparisons of differently-typed ** 2D objects. ** ***********************************************************************/ @@ -2340,75 +2303,70 @@ 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)); } 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)); + return fabs(line_calculate_point(line, pt) / + 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; + Point interpt; + LINE ln; /* * 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); + m = float8_slope(lseg->p[0].y, lseg->p[1].y, lseg->p[1].x, lseg->p[0].x); + line_construct_pm(&ln, 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) + if (lseg_interpt_line_internal(&interpt, lseg, &ln)) { /* yes, so use distance to the intersection point */ - result = point_dt(pt, ip); + result = point_dt(pt, &interpt); #ifdef GEODEBUG printf("dist_ps- distance is %f to intersection point is (%f,%f)\n", - result, ip->x, ip->y); + result, tmp.x, tmp.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; } @@ -2496,21 +2454,21 @@ dist_pb(PG_FUNCTION_ARGS) * 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_internal(NULL, lseg, line)) result = 0.0; else { result = dist_pl_internal(&lseg->p[0], line); d2 = dist_pl_internal(&lseg->p[1], line); /* XXX shouldn't we take the min not max? */ if (d2 > result) result = d2; } @@ -2649,284 +2607,258 @@ dist_ppoly_internal(Point *pt, POLYGON *poly) } /*--------------------------------------------------------------------- * 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) +/* + * Check if the line segment intersects with the line + * + * It sets the intersection point to *result, if it is not NULL. + */ +static bool +lseg_interpt_line_internal(Point *result, LSEG *lseg, LINE *line) { + Point interpt; 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", + printf("lseg_interpt_line- 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", + printf("lseg_interpt_line_- 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)) + if (line_interpt_internal(&interpt, &tmp, line)) { #ifdef GEODEBUG - printf("interpt_sl- intersection point is (%.*g %.*g)\n", DBL_DIG, p->x, DBL_DIG, p->y); + printf("lseg_interpt_line- intersection point is (%.*g %.*g)\n", DBL_DIG, p->x, DBL_DIG, p->y); #endif - if (on_ps_internal(p, lseg)) + if (on_ps_internal(&interpt, lseg)) { #ifdef GEODEBUG - printf("interpt_sl- intersection point is on segment\n"); + printf("lseg_interpt_line- intersection point is on segment\n"); #endif + if (result != NULL) + *result = interpt; + return true; } - else - p = NULL; } - return p; -} - -/* variant: just indicate if intersection point exists */ -static bool -has_interpt_sl(LSEG *lseg, LINE *line) -{ - Point *tmp; - - tmp = interpt_sl(lseg, line); - if (tmp) - return true; return false; } + /*--------------------------------------------------------------------- * close_ * Point of closest proximity between objects. *-------------------------------------------------------------------*/ /* close_pl - * The intersection point of a perpendicular of the line * through the point. */ 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? */ + point_construct(result, line->C, pt->y); + else if (FPzero(line->A)) /* horizontal? */ + point_construct(result, pt->x, line->C); + else { - 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 tmp; + + /* + * Drop a perpendicular and find the intersection point + * + * We need to invert the slope to get the perpendicular. + */ + line_construct_pm(&tmp, pt, line->B / line->A); + line_interpt_internal(result, &tmp, line); + } - /* 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 */ 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; + Point *result; double invm; int xh, yh; + LINE tmp; + + result = (Point *) palloc(sizeof(Point)); #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 /* 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 */ + *result = 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); + *result = lseg->p[yh]; /* above the lseg */ + else + /* point lines along (to left or right) of the vertical lseg. */ + point_construct(result, lseg->p[0].x, pt->y); - /* 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 */ + *result = 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); + *result = lseg->p[xh]; /* right of the lseg */ + else + /* point lines along (at top or below) the horiz. lseg. */ + point_construct(result, pt->x, lseg->p[0].y); - /* 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 + invm = float8_slope(lseg->p[0].y, lseg->p[1].y, lseg->p[1].x, lseg->p[0].x); + line_construct_pm(&tmp, &lseg->p[!yh], invm); /* lower edge of the * "band" */ - if (pt->y < (tmp->A * pt->x + tmp->C)) + 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 */ + *result = 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 + line_construct_pm(&tmp, &lseg->p[yh], invm); /* upper edge of the * "band" */ - if (pt->y > (tmp->A * pt->x + tmp->C)) + 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 */ + *result = 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); + line_construct_pm(&tmp, 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 (!lseg_interpt_line_internal(result, lseg, &tmp)) 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. */ 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; + double dist, + dist0, + dist1; - d = dist_ps_internal(&l1->p[0], l2); - dist = d; - memcpy(&point, &l1->p[0], sizeof(Point)); - - if ((d = dist_ps_internal(&l1->p[1], l2)) < dist) - { - dist = d; - memcpy(&point, &l1->p[1], sizeof(Point)); - } + dist0 = dist_ps_internal(&l1->p[0], l2); + dist1 = dist_ps_internal(&l1->p[1], l2); + dist = Min(dist0, dist1); 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), + PointPGetDatum(result), LsegPGetDatum(l2))); } - - if (dist_ps_internal(&l2->p[1], l1) < dist) + else 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), + PointPGetDatum(result), LsegPGetDatum(l2))); } - - if (result == NULL) - result = point_copy(&point); + else + { + result = (Point *) palloc(sizeof(Point)); + *result = l1->p[dist == dist0 ? 0 : 1]; + } PG_RETURN_POINT_P(result); } /* close_pb() * Closest point on or in box to specified point. */ Datum close_pb(PG_FUNCTION_ARGS) { @@ -2989,30 +2921,31 @@ close_pb(PG_FUNCTION_ARGS) 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_internal(result, lseg, line)) PG_RETURN_POINT_P(result); 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]); + *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(); } @@ -3022,30 +2955,31 @@ close_sl(PG_FUNCTION_ARGS) */ 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) + result = (Point *) palloc(sizeof(Point)); + + if (lseg_interpt_line_internal(result, lseg, line)) PG_RETURN_POINT_P(result); 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]); + *result = lseg->p[0]; else - result = point_copy(&lseg->p[1]); + *result = lseg->p[1]; PG_RETURN_POINT_P(result); } /* close_sb() * Closest point on or in box to line segment. */ Datum close_sb(PG_FUNCTION_ARGS) { @@ -3126,21 +3060,21 @@ close_lb(PG_FUNCTION_ARGS) /* on_pl - * Does the point satisfy the equation? */ 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(FPzero(line_calculate_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 */ Datum on_ps(PG_FUNCTION_ARGS) { @@ -3150,40 +3084,49 @@ on_ps(PG_FUNCTION_ARGS) PG_RETURN_BOOL(on_ps_internal(pt, lseg)); } static bool on_ps_internal(Point *pt, LSEG *lseg) { return FPeq(point_dt(pt, &lseg->p[0]) + point_dt(pt, &lseg->p[1]), point_dt(&lseg->p[0], &lseg->p[1])); } + +/* + * Check whether the point is in the box or on its border + */ 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); } + +/* + * Commutator of on_pb() + */ 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); } + /* on_ppath - * Whether a point lies within (on) a polyline. * If open, we have to (groan) check each segment. * (uses same algorithm as for point intersecting segment - tgl 1997-07-09) * If closed, we use the old O(n) ray method for point-in-polygon. * The ray is horizontal, from pt out to the right. * Each segment that crosses the ray counts as an * intersection; note that an endpoint or edge may touch * but not cross. * (we can do p-in-p in lg(n), but it takes preprocessing) @@ -3211,34 +3154,46 @@ 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)))); } + +/* + * 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. + */ 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, @@ -3250,21 +3205,21 @@ on_sb(PG_FUNCTION_ARGS) * 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_internal(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. * * Optimize for non-intersection by checking for box intersection first. @@ -3294,35 +3249,35 @@ inter_sb(PG_FUNCTION_ARGS) BoxPGetDatum(box))) || DatumGetBool(DirectFunctionCall2(on_pb, PointPGetDatum(&lseg->p[1]), BoxPGetDatum(box)))) PG_RETURN_BOOL(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)) + if (lseg_interpt_internal(NULL, &bseg, lseg)) PG_RETURN_BOOL(true); statlseg_construct(&bseg, &box->high, &point); - if (lseg_intersect_internal(&bseg, lseg)) + if (lseg_interpt_internal(NULL, &bseg, lseg)) PG_RETURN_BOOL(true); point.x = box->high.x; point.y = box->low.y; statlseg_construct(&bseg, &box->low, &point); - if (lseg_intersect_internal(&bseg, lseg)) + if (lseg_interpt_internal(NULL, &bseg, lseg)) PG_RETURN_BOOL(true); statlseg_construct(&bseg, &box->high, &point); - if (lseg_intersect_internal(&bseg, lseg)) + if (lseg_interpt_internal(NULL, &bseg, lseg)) PG_RETURN_BOOL(true); /* if we dropped through, no two segs intersected */ PG_RETURN_BOOL(false); } /* inter_lb() * Do line and box intersect? */ Datum @@ -3333,36 +3288,36 @@ inter_lb(PG_FUNCTION_ARGS) 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_internal(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_internal(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_internal(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_internal(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. @@ -3370,42 +3325,37 @@ inter_lb(PG_FUNCTION_ARGS) * make_bound_box - create the bounding box for the input polygon *------------------------------------------------------------------*/ /*--------------------------------------------------------------------- * Make the smallest bounding box for the given polygon. *---------------------------------------------------------------------*/ static void make_bound_box(POLYGON *poly) { int i; - double x1, - y1, - x2, - y2; + BOX *boundbox = &poly->boundbox; if (poly->npts > 0) { - x2 = x1 = poly->p[0].x; - y2 = y1 = poly->p[0].y; + boundbox->low.x = boundbox->high.x = poly->p[0].x; + boundbox->low.y = boundbox->high.y = 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; + if (poly->p[i].x < boundbox->low.x) + boundbox->low.x = poly->p[i].x; + if (poly->p[i].x > boundbox->high.x) + boundbox->high.x = poly->p[i].x; + if (poly->p[i].y < boundbox->low.y) + boundbox->low.y = poly->p[i].y; + if (poly->p[i].y > boundbox->high.y) + boundbox->high.y = poly->p[i].y; } - - box_fill(&(poly->boundbox), x1, x2, y1, y2); } else ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("cannot create bounding box for empty polygon"))); } /*------------------------------------------------------------------ * poly_in - read in the polygon from a string specification * @@ -3742,21 +3692,21 @@ poly_same(PG_FUNCTION_ARGS) *-----------------------------------------------------------------*/ Datum poly_overlap(PG_FUNCTION_ARGS) { POLYGON *polya = PG_GETARG_POLYGON_P(0); POLYGON *polyb = PG_GETARG_POLYGON_P(1); bool result; /* Quick check by bounding box */ result = (polya->npts > 0 && polyb->npts > 0 && - box_ov(&polya->boundbox, &polyb->boundbox)) ? true : false; + 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; @@ -3771,34 +3721,33 @@ poly_overlap(PG_FUNCTION_ARGS) { /* 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++) { sb.p[1] = polyb->p[ib]; - result = lseg_intersect_internal(&sa, &sb); + result = lseg_interpt_internal(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) { - 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); @@ -3818,27 +3767,26 @@ 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_internal(a, s->p)) { if (on_ps_internal(s->p + 1, &t)) return lseg_inside_poly(b, s->p + 1, poly, start); } - else if (POINTEQ(a, s->p + 1)) + else if (point_eq_internal(a, s->p + 1)) { if (on_ps_internal(s->p, &t)) return lseg_inside_poly(b, s->p, poly, start); } else if (on_ps_internal(s->p, &t)) { return lseg_inside_poly(b, s->p, poly, start); } else if (on_ps_internal(s->p + 1, &t)) { @@ -3861,50 +3809,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 (on_ps_internal(t.p + 1, &s)) 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)) { /* 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_internal(&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; /* @@ -4019,170 +3966,208 @@ 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); } 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_internal(result, p1, p2); PG_RETURN_POINT_P(result); } +static inline void +point_add_internal(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_internal(result, p1, p2); PG_RETURN_POINT_P(result); } +static inline void +point_sub_internal(Point *result, Point *pt1, Point *pt2) +{ + point_construct(result, + pt1->x - pt2->x, + pt1->y - pt2->y); +} + + 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_internal(result, p1, p2); PG_RETURN_POINT_P(result); } +static inline void +point_mul_internal(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_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_internal(result, p1, p2); PG_RETURN_POINT_P(result); } +static inline void +point_div_internal(Point *result, Point *pt1, Point *pt2) +{ + float8 div; + + div = (pt2->x * pt2->x) + (pt2->y * pt2->y); + point_construct(result, + ((pt1->x * pt2->x) + (pt1->y * pt2->y)) / div, + ((pt2->x * pt1->y) - (pt2->y * pt1->x)) / div); +} + /*********************************************************************** ** ** 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_internal(&result->high, &box->high, p); + point_add_internal(&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_internal(&result->high, &box->high, p); + point_sub_internal(&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_internal(&high, &box->high, p); + point_mul_internal(&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_internal(&high, &box->high, p); + point_div_internal(&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) { @@ -4278,83 +4263,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_internal(&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_internal(&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_internal(&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_internal(&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); @@ -4436,21 +4401,22 @@ poly_center(PG_FUNCTION_ARGS) 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) @@ -4468,22 +4434,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; @@ -4492,21 +4457,21 @@ poly_path(PG_FUNCTION_ARGS) /* * Never overflows: the old size fit in MaxAllocSize, and the new size is * smaller by a small constant. */ size = offsetof(PATH, p) + sizeof(path->p[0]) * poly->npts; path = (PATH *) palloc(size); SET_VARSIZE(path, size); path->npts = poly->npts; - path->closed = TRUE; + path->closed = true; /* prevent instability in unused pad bytes */ path->dummy = 0; for (i = 0; i < poly->npts; i++) { path->p[i].x = poly->p[i].x; path->p[i].y = poly->p[i].y; } PG_RETURN_PATH_P(path); @@ -4558,22 +4523,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\"", @@ -4657,22 +4621,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_internal(&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); @@ -4859,107 +4822,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_internal(&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_internal(&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_internal(&result->center, &circle->center, point); + result->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_internal(&result->center, &circle->center, point); + result->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) { @@ -4994,22 +4933,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); @@ -5372,118 +5311,55 @@ 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_internal(&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_internal(&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; + 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_internal(&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 if (ii == npts) - return TRUE; + return true; } } - return FALSE; -} - - -/*------------------------------------------------------------------------- - * Determine the hypotenuse. - * - * If required, x and y are swapped to make x the larger number. The - * traditional formula of x^2+y^2 is rearranged to factor x outside the - * sqrt. This allows computation of the hypotenuse for significantly - * larger values, and with a higher precision than when using the naive - * formula. In particular, this cannot overflow unless the final result - * would be out-of-range. - * - * sqrt( x^2 + y^2 ) = sqrt( x^2( 1 + y^2/x^2) ) - * = x * sqrt( 1 + y^2/x^2 ) - * = x * sqrt( 1 + y/x * y/x ) - * - * It is expected that this routine will eventually be replaced with the - * C99 hypot() function. - * - * This implementation conforms to IEEE Std 1003.1 and GLIBC, in that the - * case of hypot(inf,nan) results in INF, and not NAN. - *----------------------------------------------------------------------- - */ -double -pg_hypot(double x, double y) -{ - double yx; - - /* Handle INF and NaN properly */ - if (isinf(x) || isinf(y)) - return get_float8_infinity(); - - if (isnan(x) || isnan(y)) - return get_float8_nan(); - - /* Else, drop any minus signs */ - x = fabs(x); - y = fabs(y); - - /* Swap x and y if needed to make x the larger one */ - if (x < y) - { - double temp = x; - - x = y; - y = temp; - } - - /* - * If y is zero, the hypotenuse is x. This test saves a few cycles in - * such cases, but more importantly it also protects against - * divide-by-zero errors, since now x >= y. - */ - if (y == 0.0) - return x; - - /* Determine the hypotenuse */ - yx = y / x; - return x * sqrt(1.0 + (yx * yx)); + return false; } diff --git a/src/include/utils/geo_decls.h b/src/include/utils/geo_decls.h index 44c6381b85..c94e17af5f 100644 --- a/src/include/utils/geo_decls.h +++ b/src/include/utils/geo_decls.h @@ -43,21 +43,20 @@ #else #define FPzero(A) ((A) == 0) #define FPeq(A,B) ((A) == (B)) #define FPne(A,B) ((A) != (B)) #define FPlt(A,B) ((A) < (B)) #define FPle(A,B) ((A) <= (B)) #define FPgt(A,B) ((A) > (B)) #define FPge(A,B) ((A) >= (B)) #endif -#define HYPOT(A, B) pg_hypot(A, B) /*--------------------------------------------------------------------- * Point - (x,y) *-------------------------------------------------------------------*/ typedef struct { double x, y; } Point; @@ -166,21 +165,11 @@ typedef struct #define PolygonPGetDatum(X) PointerGetDatum(X) #define PG_GETARG_POLYGON_P(n) DatumGetPolygonP(PG_GETARG_DATUM(n)) #define PG_GETARG_POLYGON_P_COPY(n) DatumGetPolygonPCopy(PG_GETARG_DATUM(n)) #define PG_RETURN_POLYGON_P(x) return PolygonPGetDatum(x) #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 point 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); - -#endif /* GEO_DECLS_H */ +#endif /* GEO_DECLS_H */ diff --git a/src/test/regress/regress.c b/src/test/regress/regress.c index 0a123f2b39..4140f40b71 100644 --- a/src/test/regress/regress.c +++ b/src/test/regress/regress.c @@ -62,21 +62,23 @@ regress_dist_ptpath(PG_FUNCTION_ARGS) float8 result = 0.0; /* keep compiler quiet */ float8 tmp; int i; LSEG lseg; switch (path->npts) { case 0: PG_RETURN_NULL(); case 1: - result = point_dt(pt, &path->p[0]); + result = DatumGetFloat8(DirectFunctionCall2(point_distance, + PointPGetDatum(pt), + PointPGetDatum(&path->p[0]))); break; default: /* * the distance from a point to a path is the smallest distance * from the point to any of its constituent segments. */ Assert(path->npts > 1); for (i = 0; i < path->npts - 1; ++i) { @@ -280,22 +282,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(boxarea); Datum boxarea(PG_FUNCTION_ARGS) { BOX *box = PG_GETARG_BOX_P(0); double width, height; -- 2.13.5 (Apple Git-94)