diff options
Diffstat (limited to 'src/backend/utils/adt/geo_ops.c')
| -rw-r--r-- | src/backend/utils/adt/geo_ops.c | 1882 |
1 files changed, 862 insertions, 1020 deletions
diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c index f57380a4df2..3e7519015d2 100644 --- a/src/backend/utils/adt/geo_ops.c +++ b/src/backend/utils/adt/geo_ops.c @@ -38,24 +38,62 @@ 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); @@ -67,17 +105,6 @@ 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); /* @@ -93,6 +120,8 @@ static double dist_ppoly_internal(Point *pt, POLYGON *poly); #define RDELIM_EP ']' #define LDELIM_C '<' #define RDELIM_C '>' +#define LDELIM_L '{' +#define RDELIM_L '}' /* @@ -236,12 +265,9 @@ path_decode(char *str, bool opentype, int npts, Point *p, 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++; @@ -440,55 +466,29 @@ box_send(PG_FUNCTION_ARGS) /* 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; } @@ -505,10 +505,8 @@ 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? @@ -637,10 +635,7 @@ 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? @@ -651,10 +646,19 @@ 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); } @@ -757,7 +761,7 @@ 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)); } @@ -769,7 +773,7 @@ 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)); } @@ -787,7 +791,7 @@ box_distance(PG_FUNCTION_ARGS) 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)); } @@ -905,7 +909,7 @@ line_decode(char *s, const char *str, LINE *line) 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++; @@ -926,7 +930,7 @@ line_in(PG_FUNCTION_ARGS) s = str; while (isspace((unsigned char) *s)) s++; - if (*s == '{') + if (*s == LDELIM_L) { if (!line_decode(s + 1, str, line)) ereport(ERROR, @@ -940,12 +944,12 @@ line_in(PG_FUNCTION_ARGS) } 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); @@ -960,7 +964,8 @@ line_out(PG_FUNCTION_ARGS) 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)); } /* @@ -1003,14 +1008,12 @@ line_send(PG_FUNCTION_ARGS) * 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" */ @@ -1025,50 +1028,6 @@ line_construct_pm(Point *pt, double 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() @@ -1081,7 +1040,8 @@ line_construct_pp(PG_FUNCTION_ARGS) 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); } @@ -1096,9 +1056,7 @@ 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 @@ -1107,10 +1065,7 @@ 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 @@ -1169,6 +1124,20 @@ line_eq(PG_FUNCTION_ARGS) * 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. */ @@ -1178,16 +1147,14 @@ 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)) /* intersecting? */ 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); } @@ -1201,9 +1168,9 @@ line_interpt(PG_FUNCTION_ARGS) 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); } @@ -1211,27 +1178,25 @@ line_interpt(PG_FUNCTION_ARGS) /* * 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); } @@ -1242,18 +1207,17 @@ line_interpt_internal(LINE *l1, LINE *l2) } 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); -#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 + if (result != NULL) + point_construct(result, x, y); - return result; + return true; } @@ -1562,8 +1526,7 @@ path_inter(PG_FUNCTION_ARGS) 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; @@ -1615,7 +1578,7 @@ path_inter(PG_FUNCTION_ARGS) 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); } } @@ -1670,9 +1633,7 @@ path_distance(PG_FUNCTION_ARGS) 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; @@ -1780,30 +1741,14 @@ point_send(PG_FUNCTION_ARGS) } -static Point * -point_construct(double x, double y) +/* + * Initialize a point + */ +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; } @@ -1876,7 +1821,7 @@ 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 @@ -1885,9 +1830,16 @@ 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. *---------------------------------------------------------*/ @@ -1898,16 +1850,12 @@ 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); } @@ -1921,12 +1869,35 @@ point_slope(PG_FUNCTION_ARGS) } -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); } @@ -1952,7 +1923,7 @@ lseg_in(PG_FUNCTION_ARGS) 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); } @@ -1962,7 +1933,7 @@ 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])); } /* @@ -2012,16 +1983,13 @@ lseg_construct(PG_FUNCTION_ARGS) 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; @@ -2030,6 +1998,27 @@ statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2) 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) { @@ -2052,25 +2041,9 @@ 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) @@ -2078,39 +2051,19 @@ 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 @@ -2136,10 +2089,8 @@ 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 @@ -2148,10 +2099,8 @@ 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 @@ -2210,33 +2159,7 @@ 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)); } @@ -2254,57 +2177,38 @@ lseg_center(PG_FUNCTION_ARGS) 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; - /* - * 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; + line_construct(&tmp, &l2->p[0], lseg_sl(l2)); + if (!lseg_interpt_line(&interpt, l1, &tmp)) + return false; /* - * If the line intersection point isn't within l1 (or equivalently l2), - * there is no valid segment intersection point at all. + * If the line intersection point isn't within 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 (!lseg_contain_point(l2, &interpt)) + return false; - /* - * 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; - } + if (result != NULL) + *result = interpt; - return result; + return true; } -/* lseg_interpt - - * Find the intersection point of two segments (if any). - */ Datum lseg_interpt(PG_FUNCTION_ARGS) { @@ -2312,10 +2216,10 @@ lseg_interpt(PG_FUNCTION_ARGS) 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); } @@ -2340,15 +2244,9 @@ 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 @@ -2359,60 +2257,7 @@ 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)); } /* @@ -2429,46 +2274,34 @@ dist_ppath(PG_FUNCTION_ARGS) int i; LSEG lseg; - switch (path->npts) - { - 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); + Assert(path->npts > 0); - /* - * the distance from a point to a path is the smallest distance - * from the point to any of its constituent segments. - */ - for (i = 0; i < path->npts; i++) - { - int iprev; + /* + * 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); } @@ -2480,15 +2313,8 @@ 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)); } /* @@ -2502,12 +2328,12 @@ dist_sl(PG_FUNCTION_ARGS) 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; @@ -2524,17 +2350,8 @@ 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)); } /* @@ -2584,11 +2401,8 @@ 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 @@ -2596,11 +2410,8 @@ 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 @@ -2624,19 +2435,19 @@ dist_ppoly_internal(Point *pt, POLYGON *poly) 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 @@ -2655,49 +2466,51 @@ dist_ppoly_internal(Point *pt, POLYGON *poly) * 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 + * + * 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 +lseg_interpt_line(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", - 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; -} + /* + * First, we promote the line segment to a line, because we know how + * to find the intersection point of two lines. If they don't have + * an intersection point, we are done. + */ + line_construct(&tmp, &lseg->p[0], lseg_sl(lseg)); + if (!line_interpt_line(&interpt, &tmp, line)) + return false; -/* variant: just indicate if intersection point exists */ -static bool -has_interpt_sl(LSEG *lseg, LINE *line) -{ - Point *tmp; + /* + * Then, we check whether the intersection point is actually on the line + * segment. + */ + if (!lseg_contain_point(lseg, &interpt)) + return false; - tmp = interpt_sl(lseg, line); - if (tmp) + 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; } /*--------------------------------------------------------------------- @@ -2705,277 +2518,217 @@ has_interpt_sl(LSEG *lseg, LINE *line) * 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; + Point closept; + LINE tmp; + + /* We drop a perpendicular to find the intersection point. */ + line_construct(&tmp, point, line_invsl(line)); + retval = line_interpt_line(&closept, line, &tmp); + Assert(retval); /* XXX: We need something better. */ + + if (result != NULL) + *result = closept; + + /* Then we calculate the distance between the points. */ + return point_dt(&closept, point); +} + Datum close_pl(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); LINE *line = PG_GETARG_LINE_P(1); Point *result; - 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. */ -Datum -close_ps(PG_FUNCTION_ARGS) +static float8 +lseg_closept_point(Point *result, LSEG *lseg, Point *pt) { - Point *pt = PG_GETARG_POINT_P(0); - LSEG *lseg = PG_GETARG_LSEG_P(1); - Point *result = NULL; - LINE *tmp; - double invm; - int xh, - yh; - -#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 */ - 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); - } + Point closept; + LINE tmp; /* - * vert. and horiz. cases are down, now check if the closest point is one - * of the end points or someplace on the lseg. + * 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); - 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); - } + if (result != NULL) + *result = closept; - /* - * 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); + return point_dt(&closept, pt); +} - /* - * 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) +Datum +close_ps(PG_FUNCTION_ARGS) +{ + Point *pt = PG_GETARG_POINT_P(0); + LSEG *lseg = PG_GETARG_LSEG_P(1); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + if (isnan(lseg_closept_point(result, lseg, pt))) PG_RETURN_NULL(); -#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. +/* + * Closest point on line segment 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_lseg(PG_FUNCTION_ARGS) +static float8 +lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2) { - LSEG *l1 = PG_GETARG_LSEG_P(0); - LSEG *l2 = PG_GETARG_LSEG_P(1); - Point *result = NULL; Point point; double dist; double d; - d = dist_ps_internal(&l1->p[0], l2); + d = lseg_closept_point(NULL, l1, &l2->p[0]); dist = d; - memcpy(&point, &l1->p[0], sizeof(Point)); + if (result != NULL) + *result = l2->p[0]; - if ((d = dist_ps_internal(&l1->p[1], l2)) < dist) + d = lseg_closept_point(NULL, l1, &l2->p[1]); + if (d < dist) { dist = d; - memcpy(&point, &l1->p[1], sizeof(Point)); + if (result != NULL) + *result = l2->p[1]; } - 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 (lseg_closept_point(&point, l2, &l1->p[0]) < dist) + d = lseg_closept_point(result, l1, &point); - 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 (lseg_closept_point(&point, l2, &l1->p[1]) < dist) + d = lseg_closept_point(result, l1, &point); - if (result == NULL) - result = point_copy(&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; + + result = (Point *) palloc(sizeof(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. * @@ -2995,16 +2748,17 @@ close_sl(PG_FUNCTION_ARGS) 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 @@ -3016,93 +2770,133 @@ close_sl(PG_FUNCTION_ARGS) 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 @@ -3123,37 +2917,55 @@ close_lb(PG_FUNCTION_ARGS) * 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 @@ -3162,8 +2974,7 @@ 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 @@ -3172,8 +2983,7 @@ 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 - @@ -3217,18 +3027,33 @@ on_ppath(PG_FUNCTION_ARGS) 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 @@ -3237,12 +3062,7 @@ 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)); } /*--------------------------------------------------------------------- @@ -3256,24 +3076,28 @@ 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; @@ -3285,42 +3109,54 @@ inter_sb(PG_FUNCTION_ARGS) /* 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? */ @@ -3339,22 +3175,22 @@ inter_lb(PG_FUNCTION_ARGS) 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 */ @@ -3381,28 +3217,26 @@ make_bound_box(POLYGON *poly) 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; } /*------------------------------------------------------------------ @@ -3746,9 +3580,10 @@ poly_overlap(PG_FUNCTION_ARGS) 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 @@ -3766,7 +3601,7 @@ poly_overlap(PG_FUNCTION_ARGS) 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]; @@ -3774,10 +3609,10 @@ poly_overlap(PG_FUNCTION_ARGS) /* 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]; } @@ -3787,10 +3622,9 @@ poly_overlap(PG_FUNCTION_ARGS) 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)); } } @@ -3824,22 +3658,21 @@ touched_lseg_inside_poly(Point *a, Point *b, LSEG *s, POLYGON *poly, int start) 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); } @@ -3867,36 +3700,35 @@ lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start) 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]; @@ -3922,39 +3754,42 @@ lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start) /*----------------------------------------------------------------- * Determine if polygon A contains polygon B. *-----------------------------------------------------------------*/ -Datum -poly_contain(PG_FUNCTION_ARGS) +static bool +poly_contain_poly(POLYGON *polya, POLYGON *polyb) { - POLYGON *polya = PG_GETARG_POLYGON_P(0); - POLYGON *polyb = PG_GETARG_POLYGON_P(1); - bool result; + int i; + LSEG s; + + Assert(polya->npts > 0 && polyb->npts > 0); /* * 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; + if (!box_contain_box(&polya->boundbox, &polyb->boundbox)) + return false; - s.p[0] = polyb->p[polyb->npts - 1]; - result = true; + s.p[0] = polyb->p[polyb->npts - 1]; - 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 + for (i = 0; i < polyb->npts; i++) { - result = false; + 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; + + result = poly_contain_poly(polya, polyb); + /* * Avoid leaking memory for toasted inputs ... needed for rtree indexes */ @@ -3971,11 +3806,20 @@ poly_contain(PG_FUNCTION_ARGS) 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); } @@ -4025,8 +3869,22 @@ construct_point(PG_FUNCTION_ARGS) { float8 x = PG_GETARG_FLOAT8(0); float8 y = PG_GETARG_FLOAT8(1); + Point *result; + + result = (Point *) palloc(sizeof(Point)); + + point_construct(result, x, y); - PG_RETURN_POINT_P(point_construct(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 @@ -4038,12 +3896,20 @@ point_add(PG_FUNCTION_ARGS) 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) { @@ -4053,12 +3919,20 @@ point_sub(PG_FUNCTION_ARGS) 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) { @@ -4068,31 +3942,39 @@ point_mul(PG_FUNCTION_ARGS) 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); } -Datum -point_div(PG_FUNCTION_ARGS) + +static inline void +point_div_point(Point *result, Point *pt1, Point *pt2) { - 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); + div = (pt2->x * pt2->x) + (pt2->y * pt2->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_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; + + result = (Point *) palloc(sizeof(Point)); + + point_div_point(result, p1, p2); PG_RETURN_POINT_P(result); } @@ -4109,8 +3991,13 @@ points_box(PG_FUNCTION_ARGS) { Point *p1 = PG_GETARG_POINT_P(0); Point *p2 = PG_GETARG_POINT_P(1); + BOX *result; + + result = (BOX *) palloc(sizeof(BOX)); - PG_RETURN_BOX_P(box_construct(p1->x, p2->x, p1->y, p2->y)); + box_construct(result, p1, p2); + + PG_RETURN_BOX_P(result); } Datum @@ -4118,11 +4005,14 @@ 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 @@ -4130,11 +4020,14 @@ 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 @@ -4143,17 +4036,15 @@ 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; + + result = (BOX *) palloc(sizeof(BOX)); - high = DatumGetPointP(DirectFunctionCall2(point_mul, - PointPGetDatum(&box->high), - PointPGetDatum(p))); - low = DatumGetPointP(DirectFunctionCall2(point_mul, - PointPGetDatum(&box->low), - PointPGetDatum(p))); + point_mul_point(&high, &box->high, p); + point_mul_point(&low, &box->low, p); - result = box_construct(high->x, low->x, high->y, low->y); + box_construct(result, &high, &low); PG_RETURN_BOX_P(result); } @@ -4164,17 +4055,15 @@ 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)); + + point_div_point(&high, &box->high, p); + point_div_point(&low, &box->low, p); - result = box_construct(high->x, low->x, high->y, low->y); + box_construct(result, &high, &low); PG_RETURN_BOX_P(result); } @@ -4284,10 +4173,7 @@ path_add_pt(PG_FUNCTION_ARGS) 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); } @@ -4300,10 +4186,7 @@ path_sub_pt(PG_FUNCTION_ARGS) 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); } @@ -4316,17 +4199,10 @@ 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); } @@ -4336,17 +4212,10 @@ 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); } @@ -4421,15 +4290,15 @@ 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); } @@ -4439,10 +4308,8 @@ 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); } @@ -4474,8 +4341,7 @@ box_poly(PG_FUNCTION_ARGS) 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); } @@ -4564,8 +4430,7 @@ circle_in(PG_FUNCTION_ARGS) while (depth > 0) { - if ((*s == RDELIM) - || ((*s == RDELIM_C) && (depth == 1))) + if ((*s == RDELIM) || ((*s == RDELIM_C) && (depth == 1))) { depth--; s++; @@ -4663,8 +4528,7 @@ circle_same(PG_FUNCTION_ARGS) 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? @@ -4737,7 +4601,8 @@ 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? @@ -4748,7 +4613,8 @@ 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)); } @@ -4865,20 +4731,6 @@ circle_ge(PG_FUNCTION_ARGS) * "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. */ @@ -4889,10 +4741,10 @@ circle_add_pt(PG_FUNCTION_ARGS) 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); } @@ -4904,10 +4756,10 @@ circle_sub_pt(PG_FUNCTION_ARGS) 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); } @@ -4922,16 +4774,11 @@ 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); } @@ -4942,16 +4789,11 @@ 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); } @@ -5000,8 +4842,8 @@ circle_distance(PG_FUNCTION_ARGS) 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); @@ -5197,42 +5039,46 @@ circle_poly(PG_FUNCTION_ARGS) 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. */ -Datum -poly_circle(PG_FUNCTION_ARGS) +static void +poly_to_circle(CIRCLE *result, POLYGON *poly) { - POLYGON *poly = PG_GETARG_POLYGON_P(0); - CIRCLE *circle; int i; - if (poly->npts < 1) - ereport(ERROR, - (errcode(ERRCODE_INVALID_PARAMETER_VALUE), - errmsg("cannot convert empty polygon to circle"))); + Assert(poly->npts > 0); - circle = (CIRCLE *) palloc(sizeof(CIRCLE)); - - circle->center.x = 0; - circle->center.y = 0; - circle->radius = 0; + result->center.x = 0; + result->center.y = 0; + result->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; + 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++) - circle->radius += point_dt(&poly->p[i], &circle->center); - circle->radius /= poly->npts; + result->radius += point_dt(&poly->p[i], &result->center); + result->radius /= poly->npts; +} - PG_RETURN_CIRCLE_P(circle); +Datum +poly_circle(PG_FUNCTION_ARGS) +{ + POLYGON *poly = PG_GETARG_POLYGON_P(0); + CIRCLE *result; + + result = (CIRCLE *) palloc(sizeof(CIRCLE)); + + poly_to_circle(result, poly); + + PG_RETURN_CIRCLE_P(result); } @@ -5268,8 +5114,7 @@ point_inside(Point *p, int npts, Point *plist) 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; @@ -5378,8 +5223,7 @@ plist_same(int npts, Point *p1, Point *p2) /* 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 */ @@ -5387,8 +5231,7 @@ plist_same(int npts, Point *p1, Point *p2) { 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); @@ -5407,8 +5250,7 @@ plist_same(int npts, Point *p1, Point *p2) { 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); |