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);
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);
/*
#define RDELIM_EP ']'
#define LDELIM_C '<'
#define RDELIM_C '>'
+#define LDELIM_L '{'
+#define RDELIM_L '}'
/*
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++;
/* box_construct - fill in a new box.
*/
-static BOX *
-box_construct(double x1, double x2, double y1, double y2)
-{
- 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)
+static inline void
+box_construct(BOX *result, Point *pt1, Point *pt2)
{
- 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;
}
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?
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?
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
- PG_RETURN_BOOL(FPge(box1->high.x, box2->high.x) &&
- FPle(box1->low.x, box2->low.x) &&
- FPge(box1->high.y, box2->high.y) &&
- FPle(box1->low.y, box2->low.y));
+ PG_RETURN_BOOL(box_contain_box(box1, box2));
+}
+
+/*
+ * Check whether the box is in the box or on its border
+ */
+static bool
+box_contain_box(BOX *box1, BOX *box2)
+{
+ return FPge(box1->high.x, box2->high.x) &&
+ FPle(box1->low.x, box2->low.x) &&
+ FPge(box1->high.y, box2->high.y) &&
+ FPle(box1->low.y, box2->low.y);
}
{
BOX *box = PG_GETARG_BOX_P(0);
- PG_RETURN_FLOAT8(box->high.x - box->low.x);
+ PG_RETURN_FLOAT8(box_wd(box));
}
{
BOX *box = PG_GETARG_BOX_P(0);
- PG_RETURN_FLOAT8(box->high.y - box->low.y);
+ PG_RETURN_FLOAT8(box_ht(box));
}
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));
}
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++;
s = str;
while (isspace((unsigned char) *s))
s++;
- if (*s == '{')
+ if (*s == LDELIM_L)
{
if (!line_decode(s + 1, str, line))
ereport(ERROR,
}
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);
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));
}
/*
* Internal form: Ax+By+C=0
*---------------------------------------------------------*/
-/* line_construct_pm()
- * point-slope
+/*
+ * Fill already-allocated LINE struct from the point and the slope
*/
-static LINE *
-line_construct_pm(Point *pt, double m)
+static inline void
+line_construct(LINE *result, Point *pt, float8 m)
{
- LINE *result = (LINE *) palloc(sizeof(LINE));
-
if (m == DBL_MAX)
{
/* vertical - use "x = C" */
result->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()
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);
}
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
- PG_RETURN_BOOL(!DatumGetBool(DirectFunctionCall2(line_parallel,
- LinePGetDatum(l1),
- LinePGetDatum(l2))));
+ PG_RETURN_BOOL(line_interpt_line(NULL, l1, l2));
}
Datum
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
- if (FPzero(l1->B))
- PG_RETURN_BOOL(FPzero(l2->B));
-
- PG_RETURN_BOOL(FPeq(l2->A, l1->A * (l2->B / l1->B)));
+ PG_RETURN_BOOL(!line_interpt_line(NULL, l1, l2));
}
Datum
* Line 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.
*/
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);
}
LINE *l2 = PG_GETARG_LINE_P(1);
Point *result;
- result = line_interpt_internal(l1, l2);
+ result = (Point *) palloc(sizeof(Point));
- if (result == NULL)
+ if (!line_interpt_line(result, l1, l2))
PG_RETURN_NULL();
PG_RETURN_POINT_P(result);
}
/*
* Internal version of line_interpt
*
- * returns a NULL pointer if no intersection point
+ * This returns true if two lines intersect (they do, if they are not
+ * parallel), false if they do not. This also sets the intersection point
+ * to *result, if it is not NULL.
+ *
+ * NOTE: If the lines are identical then we will find they are parallel
+ * and report "no intersection". This is a little weird, but since
+ * there's no *unique* intersection, maybe it's appropriate behavior.
*/
-static Point *
-line_interpt_internal(LINE *l1, LINE *l2)
+static bool
+line_interpt_line(Point *result, LINE *l1, LINE *l2)
{
- Point *result;
double x,
y;
- /*
- * NOTE: if the lines are identical then we will find they are parallel
- * and report "no intersection". This is a little weird, but since
- * there's no *unique* intersection, maybe it's appropriate behavior.
- */
- if (DatumGetBool(DirectFunctionCall2(line_parallel,
- LinePGetDatum(l1),
- LinePGetDatum(l2))))
- return NULL;
-
if (FPzero(l1->B)) /* l1 vertical? */
{
+ if (FPzero(l2->B)) /* l2 vertical? */
+ return false;
+
x = l1->C;
y = (l2->A * x + l2->C);
}
}
else
{
+ if (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;
}
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;
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);
}
}
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;
}
-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;
}
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
- PG_RETURN_BOOL(FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y));
+ PG_RETURN_BOOL(point_eq_point(pt1, pt2));
}
Datum
Point *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.
*---------------------------------------------------------*/
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);
}
}
-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);
}
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);
}
{
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]));
}
/*
Point *pt2 = PG_GETARG_POINT_P(1);
LSEG *result = (LSEG *) palloc(sizeof(LSEG));
- result->p[0].x = pt1->x;
- result->p[0].y = pt1->y;
- result->p[1].x = pt2->x;
- result->p[1].y = pt2->y;
+ statlseg_construct(result, pt1, pt2);
PG_RETURN_LSEG_P(result);
}
/* like lseg_construct, but assume space already allocated */
-static void
+static inline void
statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2)
{
lseg->p[0].x = pt1->x;
lseg->p[1].y = pt2->y;
}
+
+/*
+ * Return slope of the line segment
+ */
+static inline float8
+lseg_sl(LSEG *lseg)
+{
+ return point_sl(&lseg->p[0], &lseg->p[1]);
+}
+
+
+/*
+ * Return inverse slope of the line segment
+ */
+static inline float8
+lseg_invsl(LSEG *lseg)
+{
+ return point_invsl(&lseg->p[0], &lseg->p[1]);
+}
+
+
Datum
lseg_length(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
- PG_RETURN_BOOL(lseg_intersect_internal(l1, l2));
+ PG_RETURN_BOOL(lseg_interpt_lseg(NULL, l1, l2));
}
-static bool
-lseg_intersect_internal(LSEG *l1, LSEG *l2)
-{
- LINE ln;
- Point *interpt;
- bool retval;
-
- line_construct_pts(&ln, &l2->p[0], &l2->p[1]);
- interpt = interpt_sl(l1, &ln);
-
- if (interpt != NULL && on_ps_internal(interpt, l2))
- retval = true; /* interpt on l1 and l2 */
- else
- retval = false;
- return retval;
-}
Datum
lseg_parallel(PG_FUNCTION_ARGS)
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
- PG_RETURN_BOOL(FPeq(point_sl(&l1->p[0], &l1->p[1]),
- point_sl(&l2->p[0], &l2->p[1])));
+ PG_RETURN_BOOL(FPeq(lseg_sl(l1), lseg_sl(l2)));
}
-/* lseg_perp()
+/*
* Determine if two line segments are perpendicular.
- *
- * This code did not get the correct answer for
- * '((0,0),(0,1))'::lseg ?-| '((0,0),(1,0))'::lseg
- * So, modified it to check explicitly for slope of vertical line
- * returned by point_sl() and the results seem better.
- * - thomas 1998-01-31
*/
Datum
lseg_perp(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
- double m1,
- m2;
- m1 = point_sl(&(l1->p[0]), &(l1->p[1]));
- m2 = point_sl(&(l2->p[0]), &(l2->p[1]));
-
-#ifdef GEODEBUG
- printf("lseg_perp- slopes are %g and %g\n", m1, m2);
-#endif
- if (FPzero(m1))
- PG_RETURN_BOOL(FPeq(m2, DBL_MAX));
- else if (FPzero(m2))
- PG_RETURN_BOOL(FPeq(m1, DBL_MAX));
-
- PG_RETURN_BOOL(FPeq(m1 / m2, -1.0));
+ PG_RETURN_BOOL(FPeq(lseg_sl(l1), lseg_invsl(l2)));
}
Datum
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
- PG_RETURN_BOOL(FPeq(l1->p[0].x, l2->p[0].x) &&
- FPeq(l1->p[0].y, l2->p[0].y) &&
- FPeq(l1->p[1].x, l2->p[1].x) &&
- FPeq(l1->p[1].y, l2->p[1].y));
+ PG_RETURN_BOOL(point_eq_point(&l1->p[0], &l2->p[0]) &&
+ point_eq_point(&l1->p[1], &l2->p[1]));
}
Datum
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
- PG_RETURN_BOOL(!FPeq(l1->p[0].x, l2->p[0].x) ||
- !FPeq(l1->p[0].y, l2->p[0].y) ||
- !FPeq(l1->p[1].x, l2->p[1].x) ||
- !FPeq(l1->p[1].y, l2->p[1].y));
+ PG_RETURN_BOOL(!point_eq_point(&l1->p[0], &l2->p[0]) ||
+ !point_eq_point(&l1->p[1], &l2->p[1]));
}
Datum
LSEG *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));
}
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)
{
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);
}
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
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));
}
/*
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);
}
{
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));
}
/*
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;
{
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));
}
/*
{
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
{
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
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
* 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;
}
/*---------------------------------------------------------------------
* 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.
*
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
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)
{
* 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
Point *pt = PG_GETARG_POINT_P(0);
BOX *box = PG_GETARG_BOX_P(1);
- PG_RETURN_BOOL(pt->x <= box->high.x && pt->x >= box->low.x &&
- pt->y <= box->high.y && pt->y >= box->low.y);
+ PG_RETURN_BOOL(box_contain_point(box, pt));
}
Datum
BOX *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 -
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
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));
}
/*---------------------------------------------------------------------
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;
/* 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?
*/
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 */
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;
}
/*------------------------------------------------------------------
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
sa.p[0] = polya->p[polya->npts - 1];
result = false;
- for (ia = 0; ia < polya->npts && result == false; ia++)
+ for (ia = 0; ia < polya->npts && !result; ia++)
{
/* Second point of polya's edge is a current one */
sa.p[1] = polya->p[ia];
/* Init first of polyb's edge with last point */
sb.p[0] = polyb->p[polyb->npts - 1];
- for (ib = 0; ib < polyb->npts && result == false; ib++)
+ for (ib = 0; ib < polyb->npts && !result; ib++)
{
sb.p[1] = polyb->p[ib];
- result = lseg_intersect_internal(&sa, &sb);
+ result = lseg_interpt_lseg(NULL, &sa, &sb);
sb.p[0] = sb.p[1];
}
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));
}
}
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);
}
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];
/*-----------------------------------------------------------------
* 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
*/
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);
}
{
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
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)
{
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)
{
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);
}
{
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 *box = PG_GETARG_BOX_P(0);
Point *p = PG_GETARG_POINT_P(1);
+ BOX *result;
+
+ result = (BOX *) palloc(sizeof(BOX));
- 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)));
+ point_add_point(&result->high, &box->high, p);
+ point_add_point(&result->low, &box->low, p);
+
+ PG_RETURN_BOX_P(result);
}
Datum
{
BOX *box = PG_GETARG_BOX_P(0);
Point *p = PG_GETARG_POINT_P(1);
+ BOX *result;
- PG_RETURN_BOX_P(box_construct((box->high.x - p->x),
- (box->low.x - p->x),
- (box->high.y - p->y),
- (box->low.y - p->y)));
+ result = (BOX *) palloc(sizeof(BOX));
+
+ point_sub_point(&result->high, &box->high, p);
+ point_sub_point(&result->low, &box->low, p);
+
+ PG_RETURN_BOX_P(result);
}
Datum
BOX *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));
+
+ 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);
}
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_div,
- PointPGetDatum(&box->high),
- PointPGetDatum(p)));
- low = DatumGetPointP(DirectFunctionCall2(point_div,
- PointPGetDatum(&box->low),
- PointPGetDatum(p)));
+ 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);
}
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);
}
int i;
for (i = 0; i < path->npts; i++)
- {
- path->p[i].x -= point->x;
- path->p[i].y -= point->y;
- }
+ point_sub_point(&path->p[i], &path->p[i], point);
PG_RETURN_PATH_P(path);
}
{
PATH *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);
}
{
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);
}
poly_center(PG_FUNCTION_ARGS)
{
POLYGON *poly = PG_GETARG_POLYGON_P(0);
- Datum result;
- CIRCLE *circle;
+ Point *result;
+ CIRCLE circle;
+
+ result = (Point *) palloc(sizeof(Point));
- circle = DatumGetCircleP(DirectFunctionCall1(poly_circle,
- PolygonPGetDatum(poly)));
- result = DirectFunctionCall1(circle_center,
- CirclePGetDatum(circle));
+ poly_to_circle(&circle, poly);
+ *result = circle.center;
- PG_RETURN_DATUM(result);
+ PG_RETURN_POINT_P(result);
}
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);
}
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);
}
while (depth > 0)
{
- if ((*s == RDELIM)
- || ((*s == RDELIM_C) && (depth == 1)))
+ if ((*s == RDELIM) || ((*s == RDELIM_C) && (depth == 1)))
{
depth--;
s++;
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?
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?
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));
}
* "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.
*/
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);
}
Point *point = PG_GETARG_POINT_P(1);
CIRCLE *result;
- result = circle_copy(circle);
+ result = (CIRCLE *) palloc(sizeof(CIRCLE));
- result->center.x -= point->x;
- result->center.y -= point->y;
+ point_sub_point(&result->center, &circle->center, point);
+ result->radius = circle->radius;
PG_RETURN_CIRCLE_P(result);
}
CIRCLE *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);
}
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
Point *point = PG_GETARG_POINT_P(1);
CIRCLE *result;
- Point *p;
- result = circle_copy(circle);
+ result = (CIRCLE *) palloc(sizeof(CIRCLE));
- p = DatumGetPointP(DirectFunctionCall2(point_div,
- PointPGetDatum(&circle->center),
- PointPGetDatum(point)));
- result->center.x = p->x;
- result->center.y = p->y;
- result->radius /= HYPOT(point->x, point->y);
+ point_div_point(&result->center, &circle->center, point);
+ result->radius = circle->radius / HYPOT(point->x, point->y);
PG_RETURN_CIRCLE_P(result);
}
CIRCLE *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);
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);
}
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;
/* 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 */
{
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);
{
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);
projects / postgresql / commitdiff