mirrors/postgres
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
postgres/src/backend/utils/adt/geo_ops.c
Tom Lane 8c5463a511 Increase the maximum value of extra_float_digits to 3, and have pg_dump 
use that value when the backend is new enough to allow it.  This responds
to bug report from Keh-Cheng Chu pointing out that although 2 extra digits
should be sufficient to dump and restore float8 exactly, it is possible to
need 3 extra digits for float4 values.
2009-09-11 19:17:04 +00:00

5404 lines
118 KiB
C
Raw Blame History

/*-------------------------------------------------------------------------
*
* geo_ops.c
* 2D geometric operations
*
* Portions Copyright (c) 1996-2009, PostgreSQL Global Development Group
* Portions Copyright (c) 1994, Regents of the University of California
*
*
* IDENTIFICATION
* $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.105 2009/09/11 19:17:03 tgl Exp $
*
*-------------------------------------------------------------------------
*/
#include "postgres.h"
#include <math.h>
#include <limits.h>
#include <float.h>
#include <ctype.h>
#include "libpq/pqformat.h"
#include "utils/builtins.h"
#include "utils/geo_decls.h"
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
/*
* Internal routines
*/
static int point_inside(Point *p, int npts, Point *plist);
static int lseg_crossing(double x, double y, double px, double py);
static BOX *box_construct(double x1, double x2, double y1, double y2);
static BOX *box_copy(BOX *box);
static BOX *box_fill(BOX *result, double x1, double x2, double y1, double y2);
static bool box_ov(BOX *box1, BOX *box2);
static double box_ht(BOX *box);
static double box_wd(BOX *box);
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);
static void make_bound_box(POLYGON *poly);
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 int single_decode(char *str, float8 *x, char **ss);
static int single_encode(float8 x, char *str);
static int pair_decode(char *str, float8 *x, float8 *y, char **s);
static int pair_encode(float8 x, float8 y, char *str);
static int pair_count(char *s, char delim);
static int path_decode(int opentype, int npts, char *str, int *isopen, char **ss, Point *p);
static char *path_encode(bool closed, 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);
/*
* Delimiters for input and output strings.
* LDELIM, RDELIM, and DELIM are left, right, and separator delimiters, respectively.
* LDELIM_EP, RDELIM_EP are left and right delimiters for paths with endpoints.
*/
#define LDELIM '('
#define RDELIM ')'
#define DELIM ','
#define LDELIM_EP '['
#define RDELIM_EP ']'
#define LDELIM_C '<'
#define RDELIM_C '>'
/* Maximum number of characters printed by pair_encode() */
/* ...+3+7 : 3 accounts for extra_float_digits max value */
#define P_MAXLEN (2*(DBL_DIG+3+7)+1)
/*
* Geometric data types are composed of points.
* This code tries to support a common format throughout the data types,
* to allow for more predictable usage and data type conversion.
* The fundamental unit is the point. Other units are line segments,
* open paths, boxes, closed paths, and polygons (which should be considered
* non-intersecting closed paths).
*
* Data representation is as follows:
* point: (x,y)
* line segment: [(x1,y1),(x2,y2)]
* box: (x1,y1),(x2,y2)
* open path: [(x1,y1),...,(xn,yn)]
* closed path: ((x1,y1),...,(xn,yn))
* polygon: ((x1,y1),...,(xn,yn))
*
* For boxes, the points are opposite corners with the first point at the top right.
* For closed paths and polygons, the points should be reordered to allow
* fast and correct equality comparisons.
*
* XXX perhaps points in complex shapes should be reordered internally
* to allow faster internal operations, but should keep track of input order
* and restore that order for text output - tgl 97/01/16
*/
static int
single_decode(char *str, float8 *x, char **s)
{
char *cp;
if (!PointerIsValid(str))
return FALSE;
while (isspace((unsigned char) *str))
str++;
*x = strtod(str, &cp);
#ifdef GEODEBUG
printf("single_decode- (%x) try decoding %s to %g\n", (cp - str), str, *x);
#endif
if (cp <= str)
return FALSE;
while (isspace((unsigned char) *cp))
cp++;
if (s != NULL)
*s = cp;
return TRUE;
} /* single_decode() */
static int
single_encode(float8 x, char *str)
{
int ndig = DBL_DIG + extra_float_digits;
if (ndig < 1)
ndig = 1;
sprintf(str, "%.*g", ndig, x);
return TRUE;
} /* single_encode() */
static int
pair_decode(char *str, float8 *x, float8 *y, char **s)
{
int has_delim;
char *cp;
if (!PointerIsValid(str))
return FALSE;
while (isspace((unsigned char) *str))
str++;
if ((has_delim = (*str == LDELIM)))
str++;
while (isspace((unsigned char) *str))
str++;
*x = strtod(str, &cp);
if (cp <= str)
return FALSE;
while (isspace((unsigned char) *cp))
cp++;
if (*cp++ != DELIM)
return FALSE;
while (isspace((unsigned char) *cp))
cp++;
*y = strtod(cp, &str);
if (str <= cp)
return FALSE;
while (isspace((unsigned char) *str))
str++;
if (has_delim)
{
if (*str != RDELIM)
return FALSE;
str++;
while (isspace((unsigned char) *str))
str++;
}
if (s != NULL)
*s = str;
return TRUE;
}
static int
pair_encode(float8 x, float8 y, char *str)
{
int ndig = DBL_DIG + extra_float_digits;
if (ndig < 1)
ndig = 1;
sprintf(str, "%.*g,%.*g", ndig, x, ndig, y);
return TRUE;
}
static int
path_decode(int opentype, int npts, char *str, int *isopen, char **ss, Point *p)
{
int depth = 0;
char *s,
*cp;
int i;
s = str;
while (isspace((unsigned char) *s))
s++;
if ((*isopen = (*s == LDELIM_EP)))
{
/* no open delimiter allowed? */
if (!opentype)
return FALSE;
depth++;
s++;
while (isspace((unsigned char) *s))
s++;
}
else if (*s == LDELIM)
{
cp = (s + 1);
while (isspace((unsigned char) *cp))
cp++;
if (*cp == LDELIM)
{
#ifdef NOT_USED
/* nested delimiters with only one point? */
if (npts <= 1)
return FALSE;
#endif
depth++;
s = cp;
}
else if (strrchr(s, LDELIM) == s)
{
depth++;
s = cp;
}
}
for (i = 0; i < npts; i++)
{
if (!pair_decode(s, &(p->x), &(p->y), &s))
return FALSE;
if (*s == DELIM)
s++;
p++;
}
while (depth > 0)
{
if ((*s == RDELIM)
|| ((*s == RDELIM_EP) && (*isopen) && (depth == 1)))
{
depth--;
s++;
while (isspace((unsigned char) *s))
s++;
}
else
return FALSE;
}
*ss = s;
return TRUE;
} /* path_decode() */
static char *
path_encode(bool closed, int npts, Point *pt)
{
int size = npts * (P_MAXLEN + 3) + 2;
char *result;
char *cp;
int i;
/* Check for integer overflow */
if ((size - 2) / npts != (P_MAXLEN + 3))
ereport(ERROR,
(errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
errmsg("too many points requested")));
result = palloc(size);
cp = result;
switch (closed)
{
case TRUE:
*cp++ = LDELIM;
break;
case FALSE:
*cp++ = LDELIM_EP;
break;
default:
break;
}
for (i = 0; i < npts; i++)
{
*cp++ = LDELIM;
if (!pair_encode(pt->x, pt->y, cp))
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("could not format \"path\" value")));
cp += strlen(cp);
*cp++ = RDELIM;
*cp++ = DELIM;
pt++;
}
cp--;
switch (closed)
{
case TRUE:
*cp++ = RDELIM;
break;
case FALSE:
*cp++ = RDELIM_EP;
break;
default:
break;
}
*cp = '\0';
return result;
} /* path_encode() */
/*-------------------------------------------------------------
* pair_count - count the number of points
* allow the following notation:
* '((1,2),(3,4))'
* '(1,3,2,4)'
* require an odd number of delim characters in the string
*-------------------------------------------------------------*/
static int
pair_count(char *s, char delim)
{
int ndelim = 0;
while ((s = strchr(s, delim)) != NULL)
{
ndelim++;
s++;
}
return (ndelim % 2) ? ((ndelim + 1) / 2) : -1;
}
/***********************************************************************
**
** Routines for two-dimensional boxes.
**
***********************************************************************/
/*----------------------------------------------------------
* Formatting and conversion routines.
*---------------------------------------------------------*/
/* box_in - convert a string to internal form.
*
* External format: (two corners of box)
* "(f8, f8), (f8, f8)"
* also supports the older style "(f8, f8, f8, f8)"
*/
Datum
box_in(PG_FUNCTION_ARGS)
{
char *str = PG_GETARG_CSTRING(0);
BOX *box = (BOX *) palloc(sizeof(BOX));
int isopen;
char *s;
double x,
y;
if ((!path_decode(FALSE, 2, str, &isopen, &s, &(box->high)))
|| (*s != '\0'))
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type box: \"%s\"", str)));
/* reorder corners if necessary... */
if (box->high.x < box->low.x)
{
x = box->high.x;
box->high.x = box->low.x;
box->low.x = x;
}
if (box->high.y < box->low.y)
{
y = box->high.y;
box->high.y = box->low.y;
box->low.y = y;
}
PG_RETURN_BOX_P(box);
}
/* box_out - convert a box to external form.
*/
Datum
box_out(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
PG_RETURN_CSTRING(path_encode(-1, 2, &(box->high)));
}
/*
* box_recv - converts external binary format to box
*/
Datum
box_recv(PG_FUNCTION_ARGS)
{
StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
BOX *box;
double x,
y;
box = (BOX *) palloc(sizeof(BOX));
box->high.x = pq_getmsgfloat8(buf);
box->high.y = pq_getmsgfloat8(buf);
box->low.x = pq_getmsgfloat8(buf);
box->low.y = pq_getmsgfloat8(buf);
/* reorder corners if necessary... */
if (box->high.x < box->low.x)
{
x = box->high.x;
box->high.x = box->low.x;
box->low.x = x;
}
if (box->high.y < box->low.y)
{
y = box->high.y;
box->high.y = box->low.y;
box->low.y = y;
}
PG_RETURN_BOX_P(box);
}
/*
* box_send - converts box to binary format
*/
Datum
box_send(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
StringInfoData buf;
pq_begintypsend(&buf);
pq_sendfloat8(&buf, box->high.x);
pq_sendfloat8(&buf, box->high.y);
pq_sendfloat8(&buf, box->low.x);
pq_sendfloat8(&buf, box->low.y);
PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
}
/* box_construct - fill in a new box.
*/
static 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)
{
if (x1 > x2)
{
result->high.x = x1;
result->low.x = x2;
}
else
{
result->high.x = x2;
result->low.x = x1;
}
if (y1 > y2)
{
result->high.y = y1;
result->low.y = y2;
}
else
{
result->high.y = y2;
result->low.y = y1;
}
return result;
}
/* box_copy - copy a box
*/
static BOX *
box_copy(BOX *box)
{
BOX *result = (BOX *) palloc(sizeof(BOX));
memcpy((char *) result, (char *) box, sizeof(BOX));
return result;
}
/*----------------------------------------------------------
* Relational operators for BOXes.
* <, >, <=, >=, and == are based on box area.
*---------------------------------------------------------*/
/* box_same - are two boxes identical?
*/
Datum
box_same(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPeq(box1->high.x, box2->high.x) &&
FPeq(box1->low.x, box2->low.x) &&
FPeq(box1->high.y, box2->high.y) &&
FPeq(box1->low.y, box2->low.y));
}
/* box_overlap - does box1 overlap box2?
*/
Datum
box_overlap(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(box_ov(box1, box2));
}
static bool
box_ov(BOX *box1, BOX *box2)
{
return ((FPge(box1->high.x, box2->high.x) &&
FPle(box1->low.x, box2->high.x)) ||
(FPge(box2->high.x, box1->high.x) &&
FPle(box2->low.x, box1->high.x)))
&&
((FPge(box1->high.y, box2->high.y) &&
FPle(box1->low.y, box2->high.y)) ||
(FPge(box2->high.y, box1->high.y) &&
FPle(box2->low.y, box1->high.y)));
}
/* box_left - is box1 strictly left of box2?
*/
Datum
box_left(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPlt(box1->high.x, box2->low.x));
}
/* box_overleft - is the right edge of box1 at or left of
* the right edge of box2?
*
* This is "less than or equal" for the end of a time range,
* when time ranges are stored as rectangles.
*/
Datum
box_overleft(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPle(box1->high.x, box2->high.x));
}
/* box_right - is box1 strictly right of box2?
*/
Datum
box_right(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPgt(box1->low.x, box2->high.x));
}
/* box_overright - is the left edge of box1 at or right of
* the left edge of box2?
*
* This is "greater than or equal" for time ranges, when time ranges
* are stored as rectangles.
*/
Datum
box_overright(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPge(box1->low.x, box2->low.x));
}
/* box_below - is box1 strictly below box2?
*/
Datum
box_below(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPlt(box1->high.y, box2->low.y));
}
/* box_overbelow - is the upper edge of box1 at or below
* the upper edge of box2?
*/
Datum
box_overbelow(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPle(box1->high.y, box2->high.y));
}
/* box_above - is box1 strictly above box2?
*/
Datum
box_above(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPgt(box1->low.y, box2->high.y));
}
/* box_overabove - is the lower edge of box1 at or above
* the lower edge of box2?
*/
Datum
box_overabove(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPge(box1->low.y, box2->low.y));
}
/* box_contained - is box1 contained by box2?
*/
Datum
box_contained(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(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));
}
/* box_contain - does box1 contain box2?
*/
Datum
box_contain(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPge(box1->high.x, box2->high.x) &&
FPle(box1->low.x, box2->low.x) &&
FPge(box1->high.y, box2->high.y) &&
FPle(box1->low.y, box2->low.y));
}
/* box_positionop -
* is box1 entirely {above,below} box2?
*
* box_below_eq and box_above_eq are obsolete versions that (probably
* erroneously) accept the equal-boundaries case. Since these are not
* in sync with the box_left and box_right code, they are deprecated and
* not supported in the PG 8.1 rtree operator class extension.
*/
Datum
box_below_eq(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPle(box1->high.y, box2->low.y));
}
Datum
box_above_eq(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPge(box1->low.y, box2->high.y));
}
/* box_relop - is area(box1) relop area(box2), within
* our accuracy constraint?
*/
Datum
box_lt(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPlt(box_ar(box1), box_ar(box2)));
}
Datum
box_gt(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPgt(box_ar(box1), box_ar(box2)));
}
Datum
box_eq(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPeq(box_ar(box1), box_ar(box2)));
}
Datum
box_le(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPle(box_ar(box1), box_ar(box2)));
}
Datum
box_ge(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(FPge(box_ar(box1), box_ar(box2)));
}
/*----------------------------------------------------------
* "Arithmetic" operators on boxes.
*---------------------------------------------------------*/
/* box_area - returns the area of the box.
*/
Datum
box_area(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
PG_RETURN_FLOAT8(box_ar(box));
}
/* box_width - returns the width of the box
* (horizontal magnitude).
*/
Datum
box_width(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
PG_RETURN_FLOAT8(box->high.x - box->low.x);
}
/* box_height - returns the height of the box
* (vertical magnitude).
*/
Datum
box_height(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
PG_RETURN_FLOAT8(box->high.y - box->low.y);
}
/* box_distance - returns the distance between the
* center points of two boxes.
*/
Datum
box_distance(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
Point a,
b;
box_cn(&a, box1);
box_cn(&b, box2);
PG_RETURN_FLOAT8(HYPOT(a.x - b.x, a.y - b.y));
}
/* box_center - returns the center point of the box.
*/
Datum
box_center(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
Point *result = (Point *) palloc(sizeof(Point));
box_cn(result, box);
PG_RETURN_POINT_P(result);
}
/* box_ar - returns the area of the box.
*/
static double
box_ar(BOX *box)
{
return box_wd(box) * box_ht(box);
}
/* box_cn - stores the centerpoint of the box into *center.
*/
static void
box_cn(Point *center, BOX *box)
{
center->x = (box->high.x + box->low.x) / 2.0;
center->y = (box->high.y + box->low.y) / 2.0;
}
/* box_wd - returns the width (length) of the box
* (horizontal magnitude).
*/
static double
box_wd(BOX *box)
{
return box->high.x - box->low.x;
}
/* box_ht - returns the height of the box
* (vertical magnitude).
*/
static double
box_ht(BOX *box)
{
return box->high.y - box->low.y;
}
/*----------------------------------------------------------
* Funky operations.
*---------------------------------------------------------*/
/* box_intersect -
* returns the overlapping portion of two boxes,
* or NULL if they do not intersect.
*/
Datum
box_intersect(PG_FUNCTION_ARGS)
{
BOX *box1 = PG_GETARG_BOX_P(0);
BOX *box2 = PG_GETARG_BOX_P(1);
BOX *result;
if (!box_ov(box1, box2))
PG_RETURN_NULL();
result = (BOX *) palloc(sizeof(BOX));
result->high.x = Min(box1->high.x, box2->high.x);
result->low.x = Max(box1->low.x, box2->low.x);
result->high.y = Min(box1->high.y, box2->high.y);
result->low.y = Max(box1->low.y, box2->low.y);
PG_RETURN_BOX_P(result);
}
/* box_diagonal -
* returns a line segment which happens to be the
* positive-slope diagonal of "box".
*/
Datum
box_diagonal(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
LSEG *result = (LSEG *) palloc(sizeof(LSEG));
statlseg_construct(result, &box->high, &box->low);
PG_RETURN_LSEG_P(result);
}
/***********************************************************************
**
** Routines for 2D lines.
** Lines are not intended to be used as ADTs per se,
** but their ops are useful tools for other ADT ops. Thus,
** there are few relops.
**
***********************************************************************/
Datum
line_in(PG_FUNCTION_ARGS)
{
#ifdef ENABLE_LINE_TYPE
char *str = PG_GETARG_CSTRING(0);
#endif
LINE *line;
#ifdef ENABLE_LINE_TYPE
/* when fixed, modify "not implemented", catalog/pg_type.h and SGML */
LSEG lseg;
int isopen;
char *s;
if ((!path_decode(TRUE, 2, str, &isopen, &s, &(lseg.p[0])))
|| (*s != '\0'))
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type line: \"%s\"", str)));
line = (LINE *) palloc(sizeof(LINE));
line_construct_pts(line, &lseg.p[0], &lseg.p[1]);
#else
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("type \"line\" not yet implemented")));
line = NULL;
#endif
PG_RETURN_LINE_P(line);
}
Datum
line_out(PG_FUNCTION_ARGS)
{
#ifdef ENABLE_LINE_TYPE
LINE *line = PG_GETARG_LINE_P(0);
#endif
char *result;
#ifdef ENABLE_LINE_TYPE
/* when fixed, modify "not implemented", catalog/pg_type.h and SGML */
LSEG lseg;
if (FPzero(line->B))
{ /* vertical */
/* use "x = C" */
result->A = -1;
result->B = 0;
result->C = pt1->x;
#ifdef GEODEBUG
printf("line_out- line is vertical\n");
#endif
#ifdef NOT_USED
result->m = DBL_MAX;
#endif
}
else if (FPzero(line->A))
{ /* horizontal */
/* use "x = C" */
result->A = 0;
result->B = -1;
result->C = pt1->y;
#ifdef GEODEBUG
printf("line_out- line is horizontal\n");
#endif
#ifdef NOT_USED
result->m = 0.0;
#endif
}
else
{
}
if (FPzero(line->A)) /* horizontal? */
{
}
else if (FPzero(line->B)) /* vertical? */
{
}
else
{
}
return path_encode(TRUE, 2, (Point *) &(ls->p[0]));
#else
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("type \"line\" not yet implemented")));
result = NULL;
#endif
PG_RETURN_CSTRING(result);
}
/*
* line_recv - converts external binary format to line
*/
Datum
line_recv(PG_FUNCTION_ARGS)
{
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("type \"line\" not yet implemented")));
return 0;
}
/*
* line_send - converts line to binary format
*/
Datum
line_send(PG_FUNCTION_ARGS)
{
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("type \"line\" not yet implemented")));
return 0;
}
/*----------------------------------------------------------
* Conversion routines from one line formula to internal.
* Internal form: Ax+By+C=0
*---------------------------------------------------------*/
/* line_construct_pm()
* point-slope
*/
static LINE *
line_construct_pm(Point *pt, double m)
{
LINE *result = (LINE *) palloc(sizeof(LINE));
/* use "mx - y + yinter = 0" */
result->A = m;
result->B = -1.0;
if (m == DBL_MAX)
result->C = pt->y;
else
result->C = pt->y - m * pt->x;
#ifdef NOT_USED
result->m = m;
#endif
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 NOT_USED
line->m = DBL_MAX;
#endif
#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 NOT_USED
line->m = 0.0;
#endif
#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;
#ifdef NOT_USED
line->m = line->A;
#endif
#ifdef GEODEBUG
printf("line_construct_pts- line is neither vertical nor horizontal (diffs x=%.*g, y=%.*g\n",
DBL_DIG, (pt2->x - pt1->x), DBL_DIG, (pt2->y - pt1->y));
#endif
}
}
/* line_construct_pp()
* two points
*/
Datum
line_construct_pp(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
LINE *result = (LINE *) palloc(sizeof(LINE));
line_construct_pts(result, pt1, pt2);
PG_RETURN_LINE_P(result);
}
/*----------------------------------------------------------
* Relative position routines.
*---------------------------------------------------------*/
Datum
line_intersect(PG_FUNCTION_ARGS)
{
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
PG_RETURN_BOOL(!DatumGetBool(DirectFunctionCall2(line_parallel,
LinePGetDatum(l1),
LinePGetDatum(l2))));
}
Datum
line_parallel(PG_FUNCTION_ARGS)
{
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
#ifdef NOT_USED
PG_RETURN_BOOL(FPeq(l1->m, l2->m));
#endif
if (FPzero(l1->B))
PG_RETURN_BOOL(FPzero(l2->B));
PG_RETURN_BOOL(FPeq(l2->A, l1->A * (l2->B / l1->B)));
}
Datum
line_perp(PG_FUNCTION_ARGS)
{
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
#ifdef NOT_USED
if (l1->m)
PG_RETURN_BOOL(FPeq(l2->m / l1->m, -1.0));
else if (l2->m)
PG_RETURN_BOOL(FPeq(l1->m / l2->m, -1.0));
#endif
if (FPzero(l1->A))
PG_RETURN_BOOL(FPzero(l2->B));
else if (FPzero(l1->B))
PG_RETURN_BOOL(FPzero(l2->A));
PG_RETURN_BOOL(FPeq(((l1->A * l2->B) / (l1->B * l2->A)), -1.0));
}
Datum
line_vertical(PG_FUNCTION_ARGS)
{
LINE *line = PG_GETARG_LINE_P(0);
PG_RETURN_BOOL(FPzero(line->B));
}
Datum
line_horizontal(PG_FUNCTION_ARGS)
{
LINE *line = PG_GETARG_LINE_P(0);
PG_RETURN_BOOL(FPzero(line->A));
}
Datum
line_eq(PG_FUNCTION_ARGS)
{
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
double k;
if (!FPzero(l2->A))
k = l1->A / l2->A;
else if (!FPzero(l2->B))
k = l1->B / l2->B;
else if (!FPzero(l2->C))
k = l1->C / l2->C;
else
k = 1.0;
PG_RETURN_BOOL(FPeq(l1->A, k * l2->A) &&
FPeq(l1->B, k * l2->B) &&
FPeq(l1->C, k * l2->C));
}
/*----------------------------------------------------------
* Line arithmetic routines.
*---------------------------------------------------------*/
/* line_distance()
* Distance between two lines.
*/
Datum
line_distance(PG_FUNCTION_ARGS)
{
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
float8 result;
Point *tmp;
if (!DatumGetBool(DirectFunctionCall2(line_parallel,
LinePGetDatum(l1),
LinePGetDatum(l2))))
PG_RETURN_FLOAT8(0.0);
if (FPzero(l1->B)) /* vertical? */
PG_RETURN_FLOAT8(fabs(l1->C - l2->C));
tmp = point_construct(0.0, l1->C);
result = dist_pl_internal(tmp, l2);
PG_RETURN_FLOAT8(result);
}
/* line_interpt()
* Point where two lines l1, l2 intersect (if any)
*/
Datum
line_interpt(PG_FUNCTION_ARGS)
{
LINE *l1 = PG_GETARG_LINE_P(0);
LINE *l2 = PG_GETARG_LINE_P(1);
Point *result;
result = line_interpt_internal(l1, l2);
if (result == NULL)
PG_RETURN_NULL();
PG_RETURN_POINT_P(result);
}
/*
* Internal version of line_interpt
*
* returns a NULL pointer if no intersection point
*/
static Point *
line_interpt_internal(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;
#ifdef NOT_USED
if (FPzero(l1->B)) /* l1 vertical? */
result = point_construct(l2->m * l1->C + l2->C, l1->C);
else if (FPzero(l2->B)) /* l2 vertical? */
result = point_construct(l1->m * l2->C + l1->C, l2->C);
else
{
x = (l1->C - l2->C) / (l2->A - l1->A);
result = point_construct(x, l1->m * x + l1->C);
}
#endif
if (FPzero(l1->B)) /* l1 vertical? */
{
x = l1->C;
y = (l2->A * x + l2->C);
}
else if (FPzero(l2->B)) /* l2 vertical? */
{
x = l2->C;
y = (l1->A * x + l1->C);
}
else
{
x = (l1->C - l2->C) / (l2->A - l1->A);
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
return result;
}
/***********************************************************************
**
** Routines for 2D paths (sequences of line segments, also
** called `polylines').
**
** This is not a general package for geometric paths,
** which of course include polygons; the emphasis here
** is on (for example) usefulness in wire layout.
**
***********************************************************************/
/*----------------------------------------------------------
* String to path / path to string conversion.
* External format:
* "((xcoord, ycoord),... )"
* "[(xcoord, ycoord),... ]"
* "(xcoord, ycoord),... "
* "[xcoord, ycoord,... ]"
* Also support older format:
* "(closed, npts, xcoord, ycoord,... )"
*---------------------------------------------------------*/
Datum
path_area(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P(0);
double area = 0.0;
int i,
j;
if (!path->closed)
PG_RETURN_NULL();
for (i = 0; i < path->npts; i++)
{
j = (i + 1) % path->npts;
area += path->p[i].x * path->p[j].y;
area -= path->p[i].y * path->p[j].x;
}
area *= 0.5;
PG_RETURN_FLOAT8(area < 0.0 ? -area : area);
}
Datum
path_in(PG_FUNCTION_ARGS)
{
char *str = PG_GETARG_CSTRING(0);
PATH *path;
int isopen;
char *s;
int npts;
int size;
int depth = 0;
if ((npts = pair_count(str, ',')) <= 0)
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type path: \"%s\"", str)));
s = str;
while (isspace((unsigned char) *s))
s++;
/* skip single leading paren */
if ((*s == LDELIM) && (strrchr(s, LDELIM) == s))
{
s++;
depth++;
}
size = offsetof(PATH, p[0]) +sizeof(path->p[0]) * npts;
path = (PATH *) palloc(size);
SET_VARSIZE(path, size);
path->npts = npts;
if ((!path_decode(TRUE, npts, s, &isopen, &s, &(path->p[0])))
&& (!((depth == 0) && (*s == '\0'))) && !((depth >= 1) && (*s == RDELIM)))
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type path: \"%s\"", str)));
path->closed = (!isopen);
/* prevent instability in unused pad bytes */
path->dummy = 0;
PG_RETURN_PATH_P(path);
}
Datum
path_out(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P(0);
PG_RETURN_CSTRING(path_encode(path->closed, path->npts, path->p));
}
/*
* path_recv - converts external binary format to path
*
* External representation is closed flag (a boolean byte), int32 number
* of points, and the points.
*/
Datum
path_recv(PG_FUNCTION_ARGS)
{
StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
PATH *path;
int closed;
int32 npts;
int32 i;
int size;
closed = pq_getmsgbyte(buf);
npts = pq_getmsgint(buf, sizeof(int32));
if (npts <= 0 || npts >= (int32) ((INT_MAX - offsetof(PATH, p[0])) / sizeof(Point)))
ereport(ERROR,
(errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
errmsg("invalid number of points in external \"path\" value")));
size = offsetof(PATH, p[0]) +sizeof(path->p[0]) * npts;
path = (PATH *) palloc(size);
SET_VARSIZE(path, size);
path->npts = npts;
path->closed = (closed ? 1 : 0);
for (i = 0; i < npts; i++)
{
path->p[i].x = pq_getmsgfloat8(buf);
path->p[i].y = pq_getmsgfloat8(buf);
}
PG_RETURN_PATH_P(path);
}
/*
* path_send - converts path to binary format
*/
Datum
path_send(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P(0);
StringInfoData buf;
int32 i;
pq_begintypsend(&buf);
pq_sendbyte(&buf, path->closed ? 1 : 0);
pq_sendint(&buf, path->npts, sizeof(int32));
for (i = 0; i < path->npts; i++)
{
pq_sendfloat8(&buf, path->p[i].x);
pq_sendfloat8(&buf, path->p[i].y);
}
PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
}
/*----------------------------------------------------------
* Relational operators.
* These are based on the path cardinality,
* as stupid as that sounds.
*
* Better relops and access methods coming soon.
*---------------------------------------------------------*/
Datum
path_n_lt(PG_FUNCTION_ARGS)
{
PATH *p1 = PG_GETARG_PATH_P(0);
PATH *p2 = PG_GETARG_PATH_P(1);
PG_RETURN_BOOL(p1->npts < p2->npts);
}
Datum
path_n_gt(PG_FUNCTION_ARGS)
{
PATH *p1 = PG_GETARG_PATH_P(0);
PATH *p2 = PG_GETARG_PATH_P(1);
PG_RETURN_BOOL(p1->npts > p2->npts);
}
Datum
path_n_eq(PG_FUNCTION_ARGS)
{
PATH *p1 = PG_GETARG_PATH_P(0);
PATH *p2 = PG_GETARG_PATH_P(1);
PG_RETURN_BOOL(p1->npts == p2->npts);
}
Datum
path_n_le(PG_FUNCTION_ARGS)
{
PATH *p1 = PG_GETARG_PATH_P(0);
PATH *p2 = PG_GETARG_PATH_P(1);
PG_RETURN_BOOL(p1->npts <= p2->npts);
}
Datum
path_n_ge(PG_FUNCTION_ARGS)
{
PATH *p1 = PG_GETARG_PATH_P(0);
PATH *p2 = PG_GETARG_PATH_P(1);
PG_RETURN_BOOL(p1->npts >= p2->npts);
}
/*----------------------------------------------------------
* Conversion operators.
*---------------------------------------------------------*/
Datum
path_isclosed(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P(0);
PG_RETURN_BOOL(path->closed);
}
Datum
path_isopen(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P(0);
PG_RETURN_BOOL(!path->closed);
}
Datum
path_npoints(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P(0);
PG_RETURN_INT32(path->npts);
}
Datum
path_close(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P_COPY(0);
path->closed = TRUE;
PG_RETURN_PATH_P(path);
}
Datum
path_open(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P_COPY(0);
path->closed = FALSE;
PG_RETURN_PATH_P(path);
}
/* path_inter -
* Does p1 intersect p2 at any point?
* Use bounding boxes for a quick (O(n)) check, then do a
* O(n^2) iterative edge check.
*/
Datum
path_inter(PG_FUNCTION_ARGS)
{
PATH *p1 = PG_GETARG_PATH_P(0);
PATH *p2 = PG_GETARG_PATH_P(1);
BOX b1,
b2;
int i,
j;
LSEG seg1,
seg2;
if (p1->npts <= 0 || p2->npts <= 0)
PG_RETURN_BOOL(false);
b1.high.x = b1.low.x = p1->p[0].x;
b1.high.y = b1.low.y = p1->p[0].y;
for (i = 1; i < p1->npts; i++)
{
b1.high.x = Max(p1->p[i].x, b1.high.x);
b1.high.y = Max(p1->p[i].y, b1.high.y);
b1.low.x = Min(p1->p[i].x, b1.low.x);
b1.low.y = Min(p1->p[i].y, b1.low.y);
}
b2.high.x = b2.low.x = p2->p[0].x;
b2.high.y = b2.low.y = p2->p[0].y;
for (i = 1; i < p2->npts; i++)
{
b2.high.x = Max(p2->p[i].x, b2.high.x);
b2.high.y = Max(p2->p[i].y, b2.high.y);
b2.low.x = Min(p2->p[i].x, b2.low.x);
b2.low.y = Min(p2->p[i].y, b2.low.y);
}
if (!box_ov(&b1, &b2))
PG_RETURN_BOOL(false);
/* pairwise check lseg intersections */
for (i = 0; i < p1->npts; i++)
{
int iprev;
if (i > 0)
iprev = i - 1;
else
{
if (!p1->closed)
continue;
iprev = p1->npts - 1; /* include the closure segment */
}
for (j = 0; j < p2->npts; j++)
{
int jprev;
if (j > 0)
jprev = j - 1;
else
{
if (!p2->closed)
continue;
jprev = p2->npts - 1; /* include the closure segment */
}
statlseg_construct(&seg1, &p1->p[iprev], &p1->p[i]);
statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]);
if (lseg_intersect_internal(&seg1, &seg2))
PG_RETURN_BOOL(true);
}
}
/* if we dropped through, no two segs intersected */
PG_RETURN_BOOL(false);
}
/* path_distance()
* This essentially does a cartesian product of the lsegs in the
* two paths, and finds the min distance between any two lsegs
*/
Datum
path_distance(PG_FUNCTION_ARGS)
{
PATH *p1 = PG_GETARG_PATH_P(0);
PATH *p2 = PG_GETARG_PATH_P(1);
float8 min = 0.0; /* initialize to keep compiler quiet */
bool have_min = false;
float8 tmp;
int i,
j;
LSEG seg1,
seg2;
for (i = 0; i < p1->npts; i++)
{
int iprev;
if (i > 0)
iprev = i - 1;
else
{
if (!p1->closed)
continue;
iprev = p1->npts - 1; /* include the closure segment */
}
for (j = 0; j < p2->npts; j++)
{
int jprev;
if (j > 0)
jprev = j - 1;
else
{
if (!p2->closed)
continue;
jprev = p2->npts - 1; /* include the closure segment */
}
statlseg_construct(&seg1, &p1->p[iprev], &p1->p[i]);
statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]);
tmp = DatumGetFloat8(DirectFunctionCall2(lseg_distance,
LsegPGetDatum(&seg1),
LsegPGetDatum(&seg2)));
if (!have_min || tmp < min)
{
min = tmp;
have_min = true;
}
}
}
if (!have_min)
PG_RETURN_NULL();
PG_RETURN_FLOAT8(min);
}
/*----------------------------------------------------------
* "Arithmetic" operations.
*---------------------------------------------------------*/
Datum
path_length(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P(0);
float8 result = 0.0;
int i;
for (i = 0; i < path->npts; i++)
{
int iprev;
if (i > 0)
iprev = i - 1;
else
{
if (!path->closed)
continue;
iprev = path->npts - 1; /* include the closure segment */
}
result += point_dt(&path->p[iprev], &path->p[i]);
}
PG_RETURN_FLOAT8(result);
}
/***********************************************************************
**
** Routines for 2D points.
**
***********************************************************************/
/*----------------------------------------------------------
* String to point, point to string conversion.
* External format:
* "(x,y)"
* "x,y"
*---------------------------------------------------------*/
Datum
point_in(PG_FUNCTION_ARGS)
{
char *str = PG_GETARG_CSTRING(0);
Point *point;
double x,
y;
char *s;
if (!pair_decode(str, &x, &y, &s) || (*s != '\0'))
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type point: \"%s\"", str)));
point = (Point *) palloc(sizeof(Point));
point->x = x;
point->y = y;
PG_RETURN_POINT_P(point);
}
Datum
point_out(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
PG_RETURN_CSTRING(path_encode(-1, 1, pt));
}
/*
* point_recv - converts external binary format to point
*/
Datum
point_recv(PG_FUNCTION_ARGS)
{
StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
Point *point;
point = (Point *) palloc(sizeof(Point));
point->x = pq_getmsgfloat8(buf);
point->y = pq_getmsgfloat8(buf);
PG_RETURN_POINT_P(point);
}
/*
* point_send - converts point to binary format
*/
Datum
point_send(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
StringInfoData buf;
pq_begintypsend(&buf);
pq_sendfloat8(&buf, pt->x);
pq_sendfloat8(&buf, pt->y);
PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
}
static Point *
point_construct(double x, double y)
{
Point *result = (Point *) palloc(sizeof(Point));
result->x = x;
result->y = y;
return result;
}
static Point *
point_copy(Point *pt)
{
Point *result;
if (!PointerIsValid(pt))
return NULL;
result = (Point *) palloc(sizeof(Point));
result->x = pt->x;
result->y = pt->y;
return result;
}
/*----------------------------------------------------------
* Relational operators for Points.
* Since we do have a sense of coordinates being
* "equal" to a given accuracy (point_vert, point_horiz),
* the other ops must preserve that sense. This means
* that results may, strictly speaking, be a lie (unless
* EPSILON = 0.0).
*---------------------------------------------------------*/
Datum
point_left(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
PG_RETURN_BOOL(FPlt(pt1->x, pt2->x));
}
Datum
point_right(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
PG_RETURN_BOOL(FPgt(pt1->x, pt2->x));
}
Datum
point_above(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
PG_RETURN_BOOL(FPgt(pt1->y, pt2->y));
}
Datum
point_below(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
PG_RETURN_BOOL(FPlt(pt1->y, pt2->y));
}
Datum
point_vert(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
PG_RETURN_BOOL(FPeq(pt1->x, pt2->x));
}
Datum
point_horiz(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
PG_RETURN_BOOL(FPeq(pt1->y, pt2->y));
}
Datum
point_eq(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
PG_RETURN_BOOL(FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y));
}
Datum
point_ne(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
PG_RETURN_BOOL(FPne(pt1->x, pt2->x) || FPne(pt1->y, pt2->y));
}
/*----------------------------------------------------------
* "Arithmetic" operators on points.
*---------------------------------------------------------*/
Datum
point_distance(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
PG_RETURN_FLOAT8(HYPOT(pt1->x - pt2->x, pt1->y - pt2->y));
}
double
point_dt(Point *pt1, Point *pt2)
{
#ifdef GEODEBUG
printf("point_dt- segment (%f,%f),(%f,%f) length is %f\n",
pt1->x, pt1->y, pt2->x, pt2->y, HYPOT(pt1->x - pt2->x, pt1->y - pt2->y));
#endif
return HYPOT(pt1->x - pt2->x, pt1->y - pt2->y);
}
Datum
point_slope(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
PG_RETURN_FLOAT8(point_sl(pt1, pt2));
}
double
point_sl(Point *pt1, Point *pt2)
{
return (FPeq(pt1->x, pt2->x)
? (double) DBL_MAX
: (pt1->y - pt2->y) / (pt1->x - pt2->x));
}
/***********************************************************************
**
** Routines for 2D line segments.
**
***********************************************************************/
/*----------------------------------------------------------
* String to lseg, lseg to string conversion.
* External forms: "[(x1, y1), (x2, y2)]"
* "(x1, y1), (x2, y2)"
* "x1, y1, x2, y2"
* closed form ok "((x1, y1), (x2, y2))"
* (old form) "(x1, y1, x2, y2)"
*---------------------------------------------------------*/
Datum
lseg_in(PG_FUNCTION_ARGS)
{
char *str = PG_GETARG_CSTRING(0);
LSEG *lseg;
int isopen;
char *s;
lseg = (LSEG *) palloc(sizeof(LSEG));
if ((!path_decode(TRUE, 2, str, &isopen, &s, &(lseg->p[0])))
|| (*s != '\0'))
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type lseg: \"%s\"", str)));
#ifdef NOT_USED
lseg->m = point_sl(&lseg->p[0], &lseg->p[1]);
#endif
PG_RETURN_LSEG_P(lseg);
}
Datum
lseg_out(PG_FUNCTION_ARGS)
{
LSEG *ls = PG_GETARG_LSEG_P(0);
PG_RETURN_CSTRING(path_encode(FALSE, 2, (Point *) &(ls->p[0])));
}
/*
* lseg_recv - converts external binary format to lseg
*/
Datum
lseg_recv(PG_FUNCTION_ARGS)
{
StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
LSEG *lseg;
lseg = (LSEG *) palloc(sizeof(LSEG));
lseg->p[0].x = pq_getmsgfloat8(buf);
lseg->p[0].y = pq_getmsgfloat8(buf);
lseg->p[1].x = pq_getmsgfloat8(buf);
lseg->p[1].y = pq_getmsgfloat8(buf);
#ifdef NOT_USED
lseg->m = point_sl(&lseg->p[0], &lseg->p[1]);
#endif
PG_RETURN_LSEG_P(lseg);
}
/*
* lseg_send - converts lseg to binary format
*/
Datum
lseg_send(PG_FUNCTION_ARGS)
{
LSEG *ls = PG_GETARG_LSEG_P(0);
StringInfoData buf;
pq_begintypsend(&buf);
pq_sendfloat8(&buf, ls->p[0].x);
pq_sendfloat8(&buf, ls->p[0].y);
pq_sendfloat8(&buf, ls->p[1].x);
pq_sendfloat8(&buf, ls->p[1].y);
PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
}
/* lseg_construct -
* form a LSEG from two Points.
*/
Datum
lseg_construct(PG_FUNCTION_ARGS)
{
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
LSEG *result = (LSEG *) palloc(sizeof(LSEG));
result->p[0].x = pt1->x;
result->p[0].y = pt1->y;
result->p[1].x = pt2->x;
result->p[1].y = pt2->y;
#ifdef NOT_USED
result->m = point_sl(pt1, pt2);
#endif
PG_RETURN_LSEG_P(result);
}
/* like lseg_construct, but assume space already allocated */
static void
statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2)
{
lseg->p[0].x = pt1->x;
lseg->p[0].y = pt1->y;
lseg->p[1].x = pt2->x;
lseg->p[1].y = pt2->y;
#ifdef NOT_USED
lseg->m = point_sl(pt1, pt2);
#endif
}
Datum
lseg_length(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
PG_RETURN_FLOAT8(point_dt(&lseg->p[0], &lseg->p[1]));
}
/*----------------------------------------------------------
* Relative position routines.
*---------------------------------------------------------*/
/*
** find intersection of the two lines, and see if it falls on
** both segments.
*/
Datum
lseg_intersect(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
PG_RETURN_BOOL(lseg_intersect_internal(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);
#ifdef NOT_USED
PG_RETURN_BOOL(FPeq(l1->m, l2->m));
#endif
PG_RETURN_BOOL(FPeq(point_sl(&l1->p[0], &l1->p[1]),
point_sl(&l2->p[0], &l2->p[1])));
}
/* 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));
}
Datum
lseg_vertical(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
PG_RETURN_BOOL(FPeq(lseg->p[0].x, lseg->p[1].x));
}
Datum
lseg_horizontal(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
PG_RETURN_BOOL(FPeq(lseg->p[0].y, lseg->p[1].y));
}
Datum
lseg_eq(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
PG_RETURN_BOOL(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));
}
Datum
lseg_ne(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
PG_RETURN_BOOL(!FPeq(l1->p[0].x, l2->p[0].x) ||
!FPeq(l1->p[0].y, l2->p[0].y) ||
!FPeq(l1->p[1].x, l2->p[1].x) ||
!FPeq(l1->p[1].y, l2->p[1].y));
}
Datum
lseg_lt(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
PG_RETURN_BOOL(FPlt(point_dt(&l1->p[0], &l1->p[1]),
point_dt(&l2->p[0], &l2->p[1])));
}
Datum
lseg_le(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
PG_RETURN_BOOL(FPle(point_dt(&l1->p[0], &l1->p[1]),
point_dt(&l2->p[0], &l2->p[1])));
}
Datum
lseg_gt(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
PG_RETURN_BOOL(FPgt(point_dt(&l1->p[0], &l1->p[1]),
point_dt(&l2->p[0], &l2->p[1])));
}
Datum
lseg_ge(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
PG_RETURN_BOOL(FPge(point_dt(&l1->p[0], &l1->p[1]),
point_dt(&l2->p[0], &l2->p[1])));
}
/*----------------------------------------------------------
* Line arithmetic routines.
*---------------------------------------------------------*/
/* lseg_distance -
* If two segments don't intersect, then the closest
* point will be from one of the endpoints to the other
* segment.
*/
Datum
lseg_distance(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
PG_RETURN_FLOAT8(lseg_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;
}
Datum
lseg_center(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
Point *result;
result = (Point *) palloc(sizeof(Point));
result->x = (lseg->p[0].x + lseg->p[1].x) / 2.0;
result->y = (lseg->p[0].y + lseg->p[1].y) / 2.0;
PG_RETURN_POINT_P(result);
}
static Point*
lseg_interpt_internal(LSEG *l1, LSEG *l2)
{
Point *result;
LINE tmp1,
tmp2;
/*
* Find the intersection of the appropriate lines, if any.
*/
line_construct_pts(&tmp1, &l1->p[0], &l1->p[1]);
line_construct_pts(&tmp2, &l2->p[0], &l2->p[1]);
result = line_interpt_internal(&tmp1, &tmp2);
if (!PointerIsValid(result))
return NULL;
/*
* If the line intersection point isn't within l1 (or equivalently l2),
* there is no valid segment intersection point at all.
*/
if (!on_ps_internal(result, l1) ||
!on_ps_internal(result, l2))
{
pfree(result);
return NULL;
}
/*
* If there is an intersection, then check explicitly for matching
* endpoints since there may be rounding effects with annoying lsb
* residue. - tgl 1997-07-09
*/
if ((FPeq(l1->p[0].x, l2->p[0].x) && FPeq(l1->p[0].y, l2->p[0].y)) ||
(FPeq(l1->p[0].x, l2->p[1].x) && FPeq(l1->p[0].y, l2->p[1].y)))
{
result->x = l1->p[0].x;
result->y = l1->p[0].y;
}
else if ((FPeq(l1->p[1].x, l2->p[0].x) && FPeq(l1->p[1].y, l2->p[0].y)) ||
(FPeq(l1->p[1].x, l2->p[1].x) && FPeq(l1->p[1].y, l2->p[1].y)))
{
result->x = l1->p[1].x;
result->y = l1->p[1].y;
}
return result;
}
/* lseg_interpt -
* Find the intersection point of two segments (if any).
*/
Datum
lseg_interpt(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
Point *result;
result = lseg_interpt_internal(l1, l2);
if (!PointerIsValid(result))
PG_RETURN_NULL();
PG_RETURN_POINT_P(result);
}
/***********************************************************************
**
** Routines for position comparisons of differently-typed
** 2D objects.
**
***********************************************************************/
/*---------------------------------------------------------------------
* dist_
* Minimum distance from one object to another.
*-------------------------------------------------------------------*/
Datum
dist_pl(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
LINE *line = PG_GETARG_LINE_P(1);
PG_RETURN_FLOAT8(dist_pl_internal(pt, line));
}
static double
dist_pl_internal(Point *pt, LINE *line)
{
return (line->A * pt->x + line->B * pt->y + line->C) /
HYPOT(line->A, line->B);
}
Datum
dist_ps(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
LSEG *lseg = PG_GETARG_LSEG_P(1);
PG_RETURN_FLOAT8(dist_ps_internal(pt, lseg));
}
static double
dist_ps_internal(Point *pt, LSEG *lseg)
{
double m; /* slope of perp. */
LINE *ln;
double result,
tmpdist;
Point *ip;
/*
* Construct a line perpendicular to the input segment
* and through the input point
*/
if (lseg->p[1].x == lseg->p[0].x)
m = 0;
else if (lseg->p[1].y == lseg->p[0].y)
m = (double) DBL_MAX; /* slope is infinite */
else
m = (lseg->p[0].x - lseg->p[1].x) / (lseg->p[1].y - lseg->p[0].y);
ln = line_construct_pm(pt, m);
#ifdef GEODEBUG
printf("dist_ps- line is A=%g B=%g C=%g from (point) slope (%f,%f) %g\n",
ln->A, ln->B, ln->C, pt->x, pt->y, m);
#endif
/*
* Calculate distance to the line segment or to the nearest endpoint of
* the segment.
*/
/* intersection is on the line segment? */
if ((ip = interpt_sl(lseg, ln)) != NULL)
{
/* yes, so use distance to the intersection point */
result = point_dt(pt, ip);
#ifdef GEODEBUG
printf("dist_ps- distance is %f to intersection point is (%f,%f)\n",
result, ip->x, ip->y);
#endif
}
else
{
/* no, so use distance to the nearer endpoint */
result = point_dt(pt, &lseg->p[0]);
tmpdist = point_dt(pt, &lseg->p[1]);
if (tmpdist < result)
result = tmpdist;
}
return result;
}
/*
** Distance from a point to a path
*/
Datum
dist_ppath(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
PATH *path = PG_GETARG_PATH_P(1);
float8 result = 0.0; /* keep compiler quiet */
bool have_min = false;
float8 tmp;
int i;
LSEG lseg;
switch (path->npts)
{
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);
/*
* the distance from a point to a path is the smallest distance
* from the point to any of its constituent segments.
*/
for (i = 0; i < path->npts; i++)
{
int iprev;
if (i > 0)
iprev = i - 1;
else
{
if (!path->closed)
continue;
iprev = path->npts - 1; /* include the closure segment */
}
statlseg_construct(&lseg, &path->p[iprev], &path->p[i]);
tmp = dist_ps_internal(pt, &lseg);
if (!have_min || tmp < result)
{
result = tmp;
have_min = true;
}
}
break;
}
PG_RETURN_FLOAT8(result);
}
Datum
dist_pb(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
BOX *box = PG_GETARG_BOX_P(1);
float8 result;
Point *near;
near = DatumGetPointP(DirectFunctionCall2(close_pb,
PointPGetDatum(pt),
BoxPGetDatum(box)));
result = point_dt(near, pt);
PG_RETURN_FLOAT8(result);
}
Datum
dist_sl(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
LINE *line = PG_GETARG_LINE_P(1);
float8 result,
d2;
if (has_interpt_sl(lseg, line))
result = 0.0;
else
{
result = dist_pl_internal(&lseg->p[0], line);
d2 = dist_pl_internal(&lseg->p[1], line);
/* XXX shouldn't we take the min not max? */
if (d2 > result)
result = d2;
}
PG_RETURN_FLOAT8(result);
}
Datum
dist_sb(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
BOX *box = PG_GETARG_BOX_P(1);
Point *tmp;
Datum result;
tmp = DatumGetPointP(DirectFunctionCall2(close_sb,
LsegPGetDatum(lseg),
BoxPGetDatum(box)));
result = DirectFunctionCall2(dist_pb,
PointPGetDatum(tmp),
BoxPGetDatum(box));
PG_RETURN_DATUM(result);
}
Datum
dist_lb(PG_FUNCTION_ARGS)
{
#ifdef NOT_USED
LINE *line = PG_GETARG_LINE_P(0);
BOX *box = PG_GETARG_BOX_P(1);
#endif
/* need to think about this one for a while */
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("function \"dist_lb\" not implemented")));
PG_RETURN_NULL();
}
Datum
dist_cpoly(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
POLYGON *poly = PG_GETARG_POLYGON_P(1);
float8 result;
float8 d;
int i;
LSEG seg;
if (point_inside(&(circle->center), poly->npts, poly->p) != 0)
{
#ifdef GEODEBUG
printf("dist_cpoly- center inside of polygon\n");
#endif
PG_RETURN_FLOAT8(0.0);
}
/* initialize distance with segment between first and last points */
seg.p[0].x = poly->p[0].x;
seg.p[0].y = poly->p[0].y;
seg.p[1].x = poly->p[poly->npts - 1].x;
seg.p[1].y = poly->p[poly->npts - 1].y;
result = dist_ps_internal(&circle->center, &seg);
#ifdef GEODEBUG
printf("dist_cpoly- segment 0/n distance is %f\n", result);
#endif
/* check distances for other segments */
for (i = 0; (i < poly->npts - 1); i++)
{
seg.p[0].x = poly->p[i].x;
seg.p[0].y = poly->p[i].y;
seg.p[1].x = poly->p[i + 1].x;
seg.p[1].y = poly->p[i + 1].y;
d = dist_ps_internal(&circle->center, &seg);
#ifdef GEODEBUG
printf("dist_cpoly- segment %d distance is %f\n", (i + 1), d);
#endif
if (d < result)
result = d;
}
result -= circle->radius;
if (result < 0)
result = 0;
PG_RETURN_FLOAT8(result);
}
/*---------------------------------------------------------------------
* interpt_
* Intersection point of objects.
* We choose to ignore the "point" of intersection between
* lines and boxes, since there are typically two.
*-------------------------------------------------------------------*/
/* Get intersection point of lseg and line; returns NULL if no intersection */
static Point *
interpt_sl(LSEG *lseg, LINE *line)
{
LINE tmp;
Point *p;
line_construct_pts(&tmp, &lseg->p[0], &lseg->p[1]);
p = line_interpt_internal(&tmp, line);
#ifdef GEODEBUG
printf("interpt_sl- segment is (%.*g %.*g) (%.*g %.*g)\n",
DBL_DIG, lseg->p[0].x, DBL_DIG, lseg->p[0].y, DBL_DIG, lseg->p[1].x, DBL_DIG, lseg->p[1].y);
printf("interpt_sl- segment becomes line A=%.*g B=%.*g C=%.*g\n",
DBL_DIG, tmp.A, DBL_DIG, tmp.B, DBL_DIG, tmp.C);
#endif
if (PointerIsValid(p))
{
#ifdef GEODEBUG
printf("interpt_sl- intersection point is (%.*g %.*g)\n", DBL_DIG, p->x, DBL_DIG, p->y);
#endif
if (on_ps_internal(p, lseg))
{
#ifdef GEODEBUG
printf("interpt_sl- intersection point is on segment\n");
#endif
}
else
p = NULL;
}
return p;
}
/* variant: just indicate if intersection point exists */
static bool
has_interpt_sl(LSEG *lseg, LINE *line)
{
Point *tmp;
tmp = interpt_sl(lseg, line);
if (tmp)
return true;
return false;
}
/*---------------------------------------------------------------------
* close_
* Point of closest proximity between objects.
*-------------------------------------------------------------------*/
/* close_pl -
* The intersection point of a perpendicular of the line
* through the point.
*/
Datum
close_pl(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
LINE *line = PG_GETARG_LINE_P(1);
Point *result;
LINE *tmp;
double invm;
result = (Point *) palloc(sizeof(Point));
#ifdef NOT_USED
if (FPeq(line->A, -1.0) && FPzero(line->B))
{ /* vertical */
}
#endif
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 */
#ifdef NOT_USED
invm = -1.0 / line->m;
#endif
/* invert and flip the sign on the slope to get a perpendicular */
invm = line->B / line->A;
tmp = line_construct_pm(pt, invm);
result = line_interpt_internal(tmp, line);
Assert(result != NULL);
PG_RETURN_POINT_P(result);
}
/* close_ps()
* Closest point on line segment to specified point.
* Take the closest endpoint if the point is left, right,
* above, or below the segment, otherwise find the intersection
* point of the segment and its perpendicular through the point.
*
* Some tricky code here, relying on boolean expressions
* evaluating to only zero or one to use as an array index.
* bug fixes by gthaker@atl.lmco.com; May 1, 1998
*/
Datum
close_ps(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
LSEG *lseg = PG_GETARG_LSEG_P(1);
Point *result = NULL;
LINE *tmp;
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);
}
/*
* vert. and horiz. cases are down, now check if the closest point is one
* of the end points or someplace on the lseg.
*/
invm = -1.0 / point_sl(&(lseg->p[0]), &(lseg->p[1]));
tmp = line_construct_pm(&lseg->p[!yh], invm); /* lower edge of the
* "band" */
if (pt->y < (tmp->A * pt->x + tmp->C))
{ /* we are below the lower edge */
result = point_copy(&lseg->p[!yh]); /* below the lseg, take lower
* end pt */
#ifdef GEODEBUG
printf("close_ps below: tmp A %f B %f C %f m %f\n",
tmp->A, tmp->B, tmp->C, tmp->m);
#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 m %f\n",
tmp->A, tmp->B, tmp->C, tmp->m);
#endif
PG_RETURN_POINT_P(result);
}
/*
* at this point the "normal" from point will hit lseg. The closet 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 m %f\n",
tmp->A, tmp->B, tmp->C, tmp->m);
#endif
result = interpt_sl(lseg, tmp);
Assert(result != NULL);
#ifdef GEODEBUG
printf("close_ps- result.x %f result.y %f\n", result->x, result->y);
#endif
PG_RETURN_POINT_P(result);
}
/* close_lseg()
* Closest point to l1 on l2.
*/
Datum
close_lseg(PG_FUNCTION_ARGS)
{
LSEG *l1 = PG_GETARG_LSEG_P(0);
LSEG *l2 = PG_GETARG_LSEG_P(1);
Point *result = NULL;
Point point;
double dist;
double d;
d = dist_ps_internal(&l1->p[0], l2);
dist = d;
memcpy(&point, &l1->p[0], sizeof(Point));
if ((d = dist_ps_internal(&l1->p[1], l2)) < dist)
{
dist = d;
memcpy(&point, &l1->p[1], sizeof(Point));
}
if (dist_ps_internal(&l2->p[0], l1) < dist)
{
result = DatumGetPointP(DirectFunctionCall2(close_ps,
PointPGetDatum(&l2->p[0]),
LsegPGetDatum(l1)));
memcpy(&point, result, sizeof(Point));
result = DatumGetPointP(DirectFunctionCall2(close_ps,
PointPGetDatum(&point),
LsegPGetDatum(l2)));
}
if (dist_ps_internal(&l2->p[1], l1) < dist)
{
result = DatumGetPointP(DirectFunctionCall2(close_ps,
PointPGetDatum(&l2->p[1]),
LsegPGetDatum(l1)));
memcpy(&point, result, sizeof(Point));
result = DatumGetPointP(DirectFunctionCall2(close_ps,
PointPGetDatum(&point),
LsegPGetDatum(l2)));
}
if (result == NULL)
result = point_copy(&point);
PG_RETURN_POINT_P(result);
}
/* close_pb()
* Closest point on or in box to specified point.
*/
Datum
close_pb(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
BOX *box = PG_GETARG_BOX_P(1);
LSEG lseg,
seg;
Point point;
double dist,
d;
if (DatumGetBool(DirectFunctionCall2(on_pb,
PointPGetDatum(pt),
BoxPGetDatum(box))))
PG_RETURN_POINT_P(pt);
/* 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);
statlseg_construct(&seg, &box->high, &point);
if ((d = dist_ps_internal(pt, &seg)) < dist)
{
dist = d;
memcpy(&lseg, &seg, sizeof(lseg));
}
point.x = box->high.x;
point.y = box->low.y;
statlseg_construct(&seg, &box->low, &point);
if ((d = dist_ps_internal(pt, &seg)) < dist)
{
dist = d;
memcpy(&lseg, &seg, sizeof(lseg));
}
statlseg_construct(&seg, &box->high, &point);
if ((d = dist_ps_internal(pt, &seg)) < dist)
{
dist = d;
memcpy(&lseg, &seg, sizeof(lseg));
}
PG_RETURN_DATUM(DirectFunctionCall2(close_ps,
PointPGetDatum(pt),
LsegPGetDatum(&lseg)));
}
/* close_sl()
* Closest point on line to line segment.
*
* XXX THIS CODE IS WRONG
* The code is actually calculating the point on the line segment
* which is backwards from the routine naming convention.
* Copied code to new routine close_ls() but haven't fixed this one yet.
* - thomas 1998-01-31
*/
Datum
close_sl(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
LINE *line = PG_GETARG_LINE_P(1);
Point *result;
float8 d1,
d2;
result = interpt_sl(lseg, line);
if (result)
PG_RETURN_POINT_P(result);
d1 = dist_pl_internal(&lseg->p[0], line);
d2 = dist_pl_internal(&lseg->p[1], line);
if (d1 < d2)
result = point_copy(&lseg->p[0]);
else
result = point_copy(&lseg->p[1]);
PG_RETURN_POINT_P(result);
}
/* close_ls()
* Closest point on line segment to line.
*/
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);
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]);
PG_RETURN_POINT_P(result);
}
/* close_sb()
* Closest point on or in box to line segment.
*/
Datum
close_sb(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
BOX *box = PG_GETARG_BOX_P(1);
Point point;
LSEG bseg,
seg;
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)));
}
/* 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);
statlseg_construct(&seg, &box->high, &point);
if ((d = lseg_dt(lseg, &seg)) < dist)
{
dist = d;
memcpy(&bseg, &seg, sizeof(bseg));
}
point.x = box->high.x;
point.y = box->low.y;
statlseg_construct(&seg, &box->low, &point);
if ((d = lseg_dt(lseg, &seg)) < dist)
{
dist = d;
memcpy(&bseg, &seg, sizeof(bseg));
}
statlseg_construct(&seg, &box->high, &point);
if ((d = lseg_dt(lseg, &seg)) < dist)
{
dist = d;
memcpy(&bseg, &seg, sizeof(bseg));
}
/* OK, we now have the closest line segment on the box boundary */
PG_RETURN_DATUM(DirectFunctionCall2(close_lseg,
LsegPGetDatum(lseg),
LsegPGetDatum(&bseg)));
}
Datum
close_lb(PG_FUNCTION_ARGS)
{
#ifdef NOT_USED
LINE *line = PG_GETARG_LINE_P(0);
BOX *box = PG_GETARG_BOX_P(1);
#endif
/* think about this one for a while */
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("function \"close_lb\" not implemented")));
PG_RETURN_NULL();
}
/*---------------------------------------------------------------------
* on_
* Whether one object lies completely within another.
*-------------------------------------------------------------------*/
/* on_pl -
* Does the point satisfy the equation?
*/
Datum
on_pl(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
LINE *line = PG_GETARG_LINE_P(1);
PG_RETURN_BOOL(FPzero(line->A * pt->x + line->B * pt->y + line->C));
}
/* on_ps -
* Determine colinearity by detecting a triangle inequality.
* This algorithm seems to behave nicely even with lsb residues - tgl 1997-07-09
*/
Datum
on_ps(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
LSEG *lseg = PG_GETARG_LSEG_P(1);
PG_RETURN_BOOL(on_ps_internal(pt, lseg));
}
static bool
on_ps_internal(Point *pt, LSEG *lseg)
{
return FPeq(point_dt(pt, &lseg->p[0]) + point_dt(pt, &lseg->p[1]),
point_dt(&lseg->p[0], &lseg->p[1]));
}
Datum
on_pb(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
BOX *box = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(pt->x <= box->high.x && pt->x >= box->low.x &&
pt->y <= box->high.y && pt->y >= box->low.y);
}
/* on_ppath -
* Whether a point lies within (on) a polyline.
* If open, we have to (groan) check each segment.
* (uses same algorithm as for point intersecting segment - tgl 1997-07-09)
* If closed, we use the old O(n) ray method for point-in-polygon.
* The ray is horizontal, from pt out to the right.
* Each segment that crosses the ray counts as an
* intersection; note that an endpoint or edge may touch
* but not cross.
* (we can do p-in-p in lg(n), but it takes preprocessing)
*/
Datum
on_ppath(PG_FUNCTION_ARGS)
{
Point *pt = PG_GETARG_POINT_P(0);
PATH *path = PG_GETARG_PATH_P(1);
int i,
n;
double a,
b;
/*-- OPEN --*/
if (!path->closed)
{
n = path->npts - 1;
a = point_dt(pt, &path->p[0]);
for (i = 0; i < n; i++)
{
b = point_dt(pt, &path->p[i + 1]);
if (FPeq(a + b,
point_dt(&path->p[i], &path->p[i + 1])))
PG_RETURN_BOOL(true);
a = b;
}
PG_RETURN_BOOL(false);
}
/*-- CLOSED --*/
PG_RETURN_BOOL(point_inside(pt, path->npts, path->p) != 0);
}
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))));
}
Datum
on_sb(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
BOX *box = PG_GETARG_BOX_P(1);
PG_RETURN_BOOL(DatumGetBool(DirectFunctionCall2(on_pb,
PointPGetDatum(&lseg->p[0]),
BoxPGetDatum(box))) &&
DatumGetBool(DirectFunctionCall2(on_pb,
PointPGetDatum(&lseg->p[1]),
BoxPGetDatum(box))));
}
/*---------------------------------------------------------------------
* inter_
* Whether one object intersects another.
*-------------------------------------------------------------------*/
Datum
inter_sl(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
LINE *line = PG_GETARG_LINE_P(1);
PG_RETURN_BOOL(has_interpt_sl(lseg, line));
}
/* inter_sb()
* Do line segment and box intersect?
*
* Segment completely inside box counts as intersection.
* If you want only segments crossing box boundaries,
* try converting box to path first.
*
* Optimize for non-intersection by checking for box intersection first.
* - thomas 1998-01-30
*/
Datum
inter_sb(PG_FUNCTION_ARGS)
{
LSEG *lseg = PG_GETARG_LSEG_P(0);
BOX *box = PG_GETARG_BOX_P(1);
BOX lbox;
LSEG bseg;
Point point;
lbox.low.x = Min(lseg->p[0].x, lseg->p[1].x);
lbox.low.y = Min(lseg->p[0].y, lseg->p[1].y);
lbox.high.x = Max(lseg->p[0].x, lseg->p[1].x);
lbox.high.y = Max(lseg->p[0].y, lseg->p[1].y);
/* nothing close to overlap? then not going to intersect */
if (!box_ov(&lbox, box))
PG_RETURN_BOOL(false);
/* 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);
/* 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);
statlseg_construct(&bseg, &box->high, &point);
if (lseg_intersect_internal(&bseg, lseg))
PG_RETURN_BOOL(true);
point.x = box->high.x;
point.y = box->low.y;
statlseg_construct(&bseg, &box->low, &point);
if (lseg_intersect_internal(&bseg, lseg))
PG_RETURN_BOOL(true);
statlseg_construct(&bseg, &box->high, &point);
if (lseg_intersect_internal(&bseg, lseg))
PG_RETURN_BOOL(true);
/* if we dropped through, no two segs intersected */
PG_RETURN_BOOL(false);
}
/* inter_lb()
* Do line and box intersect?
*/
Datum
inter_lb(PG_FUNCTION_ARGS)
{
LINE *line = PG_GETARG_LINE_P(0);
BOX *box = PG_GETARG_BOX_P(1);
LSEG bseg;
Point p1,
p2;
/* pairwise check lseg intersections */
p1.x = box->low.x;
p1.y = box->low.y;
p2.x = box->low.x;
p2.y = box->high.y;
statlseg_construct(&bseg, &p1, &p2);
if (has_interpt_sl(&bseg, line))
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))
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))
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))
PG_RETURN_BOOL(true);
/* if we dropped through, no intersection */
PG_RETURN_BOOL(false);
}
/*------------------------------------------------------------------
* The following routines define a data type and operator class for
* POLYGONS .... Part of which (the polygon's bounding box) is built on
* top of the BOX data type.
*
* make_bound_box - create the bounding box for the input polygon
*------------------------------------------------------------------*/
/*---------------------------------------------------------------------
* Make the smallest bounding box for the given polygon.
*---------------------------------------------------------------------*/
static void
make_bound_box(POLYGON *poly)
{
int i;
double x1,
y1,
x2,
y2;
if (poly->npts > 0)
{
x2 = x1 = poly->p[0].x;
y2 = y1 = poly->p[0].y;
for (i = 1; i < poly->npts; i++)
{
if (poly->p[i].x < x1)
x1 = poly->p[i].x;
if (poly->p[i].x > x2)
x2 = poly->p[i].x;
if (poly->p[i].y < y1)
y1 = poly->p[i].y;
if (poly->p[i].y > y2)
y2 = poly->p[i].y;
}
box_fill(&(poly->boundbox), x1, x2, y1, y2);
}
else
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("cannot create bounding box for empty polygon")));
}
/*------------------------------------------------------------------
* poly_in - read in the polygon from a string specification
*
* External format:
* "((x0,y0),...,(xn,yn))"
* "x0,y0,...,xn,yn"
* also supports the older style "(x1,...,xn,y1,...yn)"
*------------------------------------------------------------------*/
Datum
poly_in(PG_FUNCTION_ARGS)
{
char *str = PG_GETARG_CSTRING(0);
POLYGON *poly;
int npts;
int size;
int isopen;
char *s;
if ((npts = pair_count(str, ',')) <= 0)
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type polygon: \"%s\"", str)));
size = offsetof(POLYGON, p[0]) +sizeof(poly->p[0]) * npts;
poly = (POLYGON *) palloc0(size); /* zero any holes */
SET_VARSIZE(poly, size);
poly->npts = npts;
if ((!path_decode(FALSE, npts, str, &isopen, &s, &(poly->p[0])))
|| (*s != '\0'))
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type polygon: \"%s\"", str)));
make_bound_box(poly);
PG_RETURN_POLYGON_P(poly);
}
/*---------------------------------------------------------------
* poly_out - convert internal POLYGON representation to the
* character string format "((f8,f8),...,(f8,f8))"
*---------------------------------------------------------------*/
Datum
poly_out(PG_FUNCTION_ARGS)
{
POLYGON *poly = PG_GETARG_POLYGON_P(0);
PG_RETURN_CSTRING(path_encode(TRUE, poly->npts, poly->p));
}
/*
* poly_recv - converts external binary format to polygon
*
* External representation is int32 number of points, and the points.
* We recompute the bounding box on read, instead of trusting it to
* be valid. (Checking it would take just as long, so may as well
* omit it from external representation.)
*/
Datum
poly_recv(PG_FUNCTION_ARGS)
{
StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
POLYGON *poly;
int32 npts;
int32 i;
int size;
npts = pq_getmsgint(buf, sizeof(int32));
if (npts <= 0 || npts >= (int32) ((INT_MAX - offsetof(POLYGON, p[0])) / sizeof(Point)))
ereport(ERROR,
(errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
errmsg("invalid number of points in external \"polygon\" value")));
size = offsetof(POLYGON, p[0]) +sizeof(poly->p[0]) * npts;
poly = (POLYGON *) palloc0(size); /* zero any holes */
SET_VARSIZE(poly, size);
poly->npts = npts;
for (i = 0; i < npts; i++)
{
poly->p[i].x = pq_getmsgfloat8(buf);
poly->p[i].y = pq_getmsgfloat8(buf);
}
make_bound_box(poly);
PG_RETURN_POLYGON_P(poly);
}
/*
* poly_send - converts polygon to binary format
*/
Datum
poly_send(PG_FUNCTION_ARGS)
{
POLYGON *poly = PG_GETARG_POLYGON_P(0);
StringInfoData buf;
int32 i;
pq_begintypsend(&buf);
pq_sendint(&buf, poly->npts, sizeof(int32));
for (i = 0; i < poly->npts; i++)
{
pq_sendfloat8(&buf, poly->p[i].x);
pq_sendfloat8(&buf, poly->p[i].y);
}
PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
}
/*-------------------------------------------------------
* Is polygon A strictly left of polygon B? i.e. is
* the right most point of A left of the left most point
* of B?
*-------------------------------------------------------*/
Datum
poly_left(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
result = polya->boundbox.high.x < polyb->boundbox.low.x;
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
/*-------------------------------------------------------
* Is polygon A overlapping or left of polygon B? i.e. is
* the right most point of A at or left of the right most point
* of B?
*-------------------------------------------------------*/
Datum
poly_overleft(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
result = polya->boundbox.high.x <= polyb->boundbox.high.x;
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
/*-------------------------------------------------------
* Is polygon A strictly right of polygon B? i.e. is
* the left most point of A right of the right most point
* of B?
*-------------------------------------------------------*/
Datum
poly_right(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
result = polya->boundbox.low.x > polyb->boundbox.high.x;
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
/*-------------------------------------------------------
* Is polygon A overlapping or right of polygon B? i.e. is
* the left most point of A at or right of the left most point
* of B?
*-------------------------------------------------------*/
Datum
poly_overright(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
result = polya->boundbox.low.x >= polyb->boundbox.low.x;
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
/*-------------------------------------------------------
* Is polygon A strictly below polygon B? i.e. is
* the upper most point of A below the lower most point
* of B?
*-------------------------------------------------------*/
Datum
poly_below(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
result = polya->boundbox.high.y < polyb->boundbox.low.y;
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
/*-------------------------------------------------------
* Is polygon A overlapping or below polygon B? i.e. is
* the upper most point of A at or below the upper most point
* of B?
*-------------------------------------------------------*/
Datum
poly_overbelow(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
result = polya->boundbox.high.y <= polyb->boundbox.high.y;
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
/*-------------------------------------------------------
* Is polygon A strictly above polygon B? i.e. is
* the lower most point of A above the upper most point
* of B?
*-------------------------------------------------------*/
Datum
poly_above(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
result = polya->boundbox.low.y > polyb->boundbox.high.y;
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
/*-------------------------------------------------------
* Is polygon A overlapping or above polygon B? i.e. is
* the lower most point of A at or above the lower most point
* of B?
*-------------------------------------------------------*/
Datum
poly_overabove(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
result = polya->boundbox.low.y >= polyb->boundbox.low.y;
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
/*-------------------------------------------------------
* Is polygon A the same as polygon B? i.e. are all the
* points the same?
* Check all points for matches in both forward and reverse
* direction since polygons are non-directional and are
* closed shapes.
*-------------------------------------------------------*/
Datum
poly_same(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
if (polya->npts != polyb->npts)
result = false;
else
result = plist_same(polya->npts, polya->p, polyb->p);
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
/*-----------------------------------------------------------------
* Determine if polygon A overlaps polygon B
*-----------------------------------------------------------------*/
Datum
poly_overlap(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
/* Quick check by bounding box */
result = (polya->npts > 0 && polyb->npts > 0 &&
box_ov(&polya->boundbox, &polyb->boundbox)) ? true : false;
/*
* Brute-force algorithm - try to find intersected edges,
* if so then polygons are overlapped else check is one
* polygon inside other or not by testing single point
* of them.
*/
if (result)
{
int ia, ib;
LSEG sa, sb;
/* Init first of polya's edge with last point */
sa.p[0] = polya->p[polya->npts - 1];
result = false;
for(ia=0; ia<polya->npts && result == false; ia++)
{
/* Second point of polya's edge is a current one */
sa.p[1] = polya->p[ia];
/* Init first of polyb's edge with last point */
sb.p[0] = polyb->p[polyb->npts - 1];
for(ib=0; ib<polyb->npts && result == false; ib++)
{
sb.p[1] = polyb->p[ib];
result = lseg_intersect_internal(&sa, &sb);
sb.p[0] = sb.p[1];
}
/*
* move current endpoint to the first point
* of next edge
*/
sa.p[0] = sa.p[1];
}
if (result==false)
{
result = ( point_inside(polya->p, polyb->npts, polyb->p)
||
point_inside(polyb->p, polya->npts, polya->p) );
}
}
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
/*
* Tests special kind of segment for in/out of polygon.
* Special kind means:
* - point a should be on segment s
* - segment (a,b) should not be contained by s
* Returns true if:
* - segment (a,b) is collinear to s and (a,b) is in polygon
* - segment (a,b) s not collinear to s. Note: that doesn't
* mean that segment is in polygon!
*/
static bool
touched_lseg_inside_poly(Point *a, Point *b, LSEG *s, POLYGON *poly, int start)
{
/* point a is on s, b is not */
LSEG t;
t.p[0] = *a;
t.p[1] = *b;
#define POINTEQ(pt1, pt2) (FPeq((pt1)->x, (pt2)->x) && FPeq((pt1)->y, (pt2)->y))
if ( POINTEQ(a, s->p) )
{
if ( on_ps_internal(s->p+1, &t) )
return lseg_inside_poly(b, s->p+1, poly, start);
}
else if (POINTEQ(a, s->p+1))
{
if ( on_ps_internal(s->p, &t) )
return lseg_inside_poly(b, s->p, poly, start);
}
else if ( on_ps_internal(s->p, &t) )
{
return lseg_inside_poly(b, s->p, poly, start);
}
else if ( on_ps_internal(s->p+1, &t) )
{
return lseg_inside_poly(b, s->p+1, poly, start);
}
return true; /* may be not true, but that will check later */
}
/*
* Returns true if segment (a,b) is in polygon, option
* start is used for optimization - function checks
* polygon's edges started from start
*/
static bool
lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start)
{
LSEG s,
t;
int i;
bool res = true,
intersection = false;
t.p[0] = *a;
t.p[1] = *b;
s.p[0] = poly->p[( start == 0) ? (poly->npts - 1) : (start - 1)];
for(i=start; i<poly->npts && res == true; i++)
{
Point *interpt;
s.p[1] = poly->p[i];
if ( on_ps_internal(t.p, &s) )
{
if ( on_ps_internal(t.p+1, &s) )
return true; /* t is contained by s */
/* Y-cross */
res = touched_lseg_inside_poly(t.p, t.p+1, &s, poly, i+1);
}
else if ( on_ps_internal(t.p+1, &s) )
{
/* Y-cross */
res = touched_lseg_inside_poly(t.p+1, t.p, &s, poly, i+1);
}
else if ( (interpt = lseg_interpt_internal(&t, &s)) != NULL )
{
/*
* segments are X-crossing, go to check each subsegment
*/
intersection = true;
res = lseg_inside_poly(t.p, interpt, poly, i+1);
if (res)
res = lseg_inside_poly(t.p+1, interpt, poly, i+1);
pfree(interpt);
}
s.p[0] = s.p[1];
}
if (res && !intersection)
{
Point p;
/*
* if X-intersection wasn't found then check central point
* of tested segment. In opposite case we already check all
* subsegments
*/
p.x = (t.p[0].x + t.p[1].x) / 2.0;
p.y = (t.p[0].y + t.p[1].y) / 2.0;
res = point_inside(&p, poly->npts, poly->p);
}
return res;
}
/*-----------------------------------------------------------------
* Determine if polygon A contains polygon B.
*-----------------------------------------------------------------*/
Datum
poly_contain(PG_FUNCTION_ARGS)
{
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
bool result;
/*
* Quick check to see if bounding box is contained.
*/
if (polya->npts > 0 && polyb->npts > 0 &&
DatumGetBool(DirectFunctionCall2(box_contain,
BoxPGetDatum(&polya->boundbox),
BoxPGetDatum(&polyb->boundbox))))
{
int i;
LSEG s;
s.p[0] = polyb->p[polyb->npts - 1];
result = true;
for(i=0; i<polyb->npts && result == true; i++)
{
s.p[1] = polyb->p[i];
result = lseg_inside_poly(s.p, s.p+1, polya, 0);
s.p[0] = s.p[1];
}
}
else
{
result = false;
}
/*
* Avoid leaking memory for toasted inputs ... needed for rtree indexes
*/
PG_FREE_IF_COPY(polya, 0);
PG_FREE_IF_COPY(polyb, 1);
PG_RETURN_BOOL(result);
}
/*-----------------------------------------------------------------
* Determine if polygon A is contained by polygon B
*-----------------------------------------------------------------*/
Datum
poly_contained(PG_FUNCTION_ARGS)
{
Datum polya = PG_GETARG_DATUM(0);
Datum polyb = PG_GETARG_DATUM(1);
/* Just switch the arguments and pass it off to poly_contain */
PG_RETURN_DATUM(DirectFunctionCall2(poly_contain, polyb, polya));
}
Datum
poly_contain_pt(PG_FUNCTION_ARGS)
{
POLYGON *poly = PG_GETARG_POLYGON_P(0);
Point *p = PG_GETARG_POINT_P(1);
PG_RETURN_BOOL(point_inside(p, poly->npts, poly->p) != 0);
}
Datum
pt_contained_poly(PG_FUNCTION_ARGS)
{
Point *p = PG_GETARG_POINT_P(0);
POLYGON *poly = PG_GETARG_POLYGON_P(1);
PG_RETURN_BOOL(point_inside(p, poly->npts, poly->p) != 0);
}
Datum
poly_distance(PG_FUNCTION_ARGS)
{
#ifdef NOT_USED
POLYGON *polya = PG_GETARG_POLYGON_P(0);
POLYGON *polyb = PG_GETARG_POLYGON_P(1);
#endif
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("function \"poly_distance\" not implemented")));
PG_RETURN_NULL();
}
/***********************************************************************
**
** Routines for 2D points.
**
***********************************************************************/
Datum
construct_point(PG_FUNCTION_ARGS)
{
float8 x = PG_GETARG_FLOAT8(0);
float8 y = PG_GETARG_FLOAT8(1);
PG_RETURN_POINT_P(point_construct(x, y));
}
Datum
point_add(PG_FUNCTION_ARGS)
{
Point *p1 = PG_GETARG_POINT_P(0);
Point *p2 = PG_GETARG_POINT_P(1);
Point *result;
result = (Point *) palloc(sizeof(Point));
result->x = (p1->x + p2->x);
result->y = (p1->y + p2->y);
PG_RETURN_POINT_P(result);
}
Datum
point_sub(PG_FUNCTION_ARGS)
{
Point *p1 = PG_GETARG_POINT_P(0);
Point *p2 = PG_GETARG_POINT_P(1);
Point *result;
result = (Point *) palloc(sizeof(Point));
result->x = (p1->x - p2->x);
result->y = (p1->y - p2->y);
PG_RETURN_POINT_P(result);
}
Datum
point_mul(PG_FUNCTION_ARGS)
{
Point *p1 = PG_GETARG_POINT_P(0);
Point *p2 = PG_GETARG_POINT_P(1);
Point *result;
result = (Point *) palloc(sizeof(Point));
result->x = (p1->x * p2->x) - (p1->y * p2->y);
result->y = (p1->x * p2->y) + (p1->y * p2->x);
PG_RETURN_POINT_P(result);
}
Datum
point_div(PG_FUNCTION_ARGS)
{
Point *p1 = PG_GETARG_POINT_P(0);
Point *p2 = PG_GETARG_POINT_P(1);
Point *result;
double div;
result = (Point *) palloc(sizeof(Point));
div = (p2->x * p2->x) + (p2->y * p2->y);
if (div == 0.0)
ereport(ERROR,
(errcode(ERRCODE_DIVISION_BY_ZERO),
errmsg("division by zero")));
result->x = ((p1->x * p2->x) + (p1->y * p2->y)) / div;
result->y = ((p2->x * p1->y) - (p2->y * p1->x)) / div;
PG_RETURN_POINT_P(result);
}
/***********************************************************************
**
** Routines for 2D boxes.
**
***********************************************************************/
Datum
points_box(PG_FUNCTION_ARGS)
{
Point *p1 = PG_GETARG_POINT_P(0);
Point *p2 = PG_GETARG_POINT_P(1);
PG_RETURN_BOX_P(box_construct(p1->x, p2->x, p1->y, p2->y));
}
Datum
box_add(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
Point *p = PG_GETARG_POINT_P(1);
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)));
}
Datum
box_sub(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
Point *p = PG_GETARG_POINT_P(1);
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)));
}
Datum
box_mul(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
Point *p = PG_GETARG_POINT_P(1);
BOX *result;
Point *high,
*low;
high = DatumGetPointP(DirectFunctionCall2(point_mul,
PointPGetDatum(&box->high),
PointPGetDatum(p)));
low = DatumGetPointP(DirectFunctionCall2(point_mul,
PointPGetDatum(&box->low),
PointPGetDatum(p)));
result = box_construct(high->x, low->x, high->y, low->y);
PG_RETURN_BOX_P(result);
}
Datum
box_div(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
Point *p = PG_GETARG_POINT_P(1);
BOX *result;
Point *high,
*low;
high = DatumGetPointP(DirectFunctionCall2(point_div,
PointPGetDatum(&box->high),
PointPGetDatum(p)));
low = DatumGetPointP(DirectFunctionCall2(point_div,
PointPGetDatum(&box->low),
PointPGetDatum(p)));
result = box_construct(high->x, low->x, high->y, low->y);
PG_RETURN_BOX_P(result);
}
/***********************************************************************
**
** Routines for 2D paths.
**
***********************************************************************/
/* path_add()
* Concatenate two paths (only if they are both open).
*/
Datum
path_add(PG_FUNCTION_ARGS)
{
PATH *p1 = PG_GETARG_PATH_P(0);
PATH *p2 = PG_GETARG_PATH_P(1);
PATH *result;
int size,
base_size;
int i;
if (p1->closed || p2->closed)
PG_RETURN_NULL();
base_size = sizeof(p1->p[0]) * (p1->npts + p2->npts);
size = offsetof(PATH, p[0]) +base_size;
/* Check for integer overflow */
if (base_size / sizeof(p1->p[0]) != (p1->npts + p2->npts) ||
size <= base_size)
ereport(ERROR,
(errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
errmsg("too many points requested")));
result = (PATH *) palloc(size);
SET_VARSIZE(result, size);
result->npts = (p1->npts + p2->npts);
result->closed = p1->closed;
for (i = 0; i < p1->npts; i++)
{
result->p[i].x = p1->p[i].x;
result->p[i].y = p1->p[i].y;
}
for (i = 0; i < p2->npts; i++)
{
result->p[i + p1->npts].x = p2->p[i].x;
result->p[i + p1->npts].y = p2->p[i].y;
}
PG_RETURN_PATH_P(result);
}
/* path_add_pt()
* Translation operators.
*/
Datum
path_add_pt(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P_COPY(0);
Point *point = PG_GETARG_POINT_P(1);
int i;
for (i = 0; i < path->npts; i++)
{
path->p[i].x += point->x;
path->p[i].y += point->y;
}
PG_RETURN_PATH_P(path);
}
Datum
path_sub_pt(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P_COPY(0);
Point *point = PG_GETARG_POINT_P(1);
int i;
for (i = 0; i < path->npts; i++)
{
path->p[i].x -= point->x;
path->p[i].y -= point->y;
}
PG_RETURN_PATH_P(path);
}
/* path_mul_pt()
* Rotation and scaling operators.
*/
Datum
path_mul_pt(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P_COPY(0);
Point *point = PG_GETARG_POINT_P(1);
Point *p;
int i;
for (i = 0; i < path->npts; i++)
{
p = DatumGetPointP(DirectFunctionCall2(point_mul,
PointPGetDatum(&path->p[i]),
PointPGetDatum(point)));
path->p[i].x = p->x;
path->p[i].y = p->y;
}
PG_RETURN_PATH_P(path);
}
Datum
path_div_pt(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P_COPY(0);
Point *point = PG_GETARG_POINT_P(1);
Point *p;
int i;
for (i = 0; i < path->npts; i++)
{
p = DatumGetPointP(DirectFunctionCall2(point_div,
PointPGetDatum(&path->p[i]),
PointPGetDatum(point)));
path->p[i].x = p->x;
path->p[i].y = p->y;
}
PG_RETURN_PATH_P(path);
}
Datum
path_center(PG_FUNCTION_ARGS)
{
#ifdef NOT_USED
PATH *path = PG_GETARG_PATH_P(0);
#endif
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("function \"path_center\" not implemented")));
PG_RETURN_NULL();
}
Datum
path_poly(PG_FUNCTION_ARGS)
{
PATH *path = PG_GETARG_PATH_P(0);
POLYGON *poly;
int size;
int i;
/* This is not very consistent --- other similar cases return NULL ... */
if (!path->closed)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("open path cannot be converted to polygon")));
size = offsetof(POLYGON, p[0]) +sizeof(poly->p[0]) * path->npts;
poly = (POLYGON *) palloc(size);
SET_VARSIZE(poly, size);
poly->npts = path->npts;
for (i = 0; i < path->npts; i++)
{
poly->p[i].x = path->p[i].x;
poly->p[i].y = path->p[i].y;
}
make_bound_box(poly);
PG_RETURN_POLYGON_P(poly);
}
/***********************************************************************
**
** Routines for 2D polygons.
**
***********************************************************************/
Datum
poly_npoints(PG_FUNCTION_ARGS)
{
POLYGON *poly = PG_GETARG_POLYGON_P(0);
PG_RETURN_INT32(poly->npts);
}
Datum
poly_center(PG_FUNCTION_ARGS)
{
POLYGON *poly = PG_GETARG_POLYGON_P(0);
Datum result;
CIRCLE *circle;
circle = DatumGetCircleP(DirectFunctionCall1(poly_circle,
PolygonPGetDatum(poly)));
result = DirectFunctionCall1(circle_center,
CirclePGetDatum(circle));
PG_RETURN_DATUM(result);
}
Datum
poly_box(PG_FUNCTION_ARGS)
{
POLYGON *poly = PG_GETARG_POLYGON_P(0);
BOX *box;
if (poly->npts < 1)
PG_RETURN_NULL();
box = box_copy(&poly->boundbox);
PG_RETURN_BOX_P(box);
}
/* box_poly()
* Convert a box to a polygon.
*/
Datum
box_poly(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
POLYGON *poly;
int size;
/* map four corners of the box to a polygon */
size = offsetof(POLYGON, p[0]) +sizeof(poly->p[0]) * 4;
poly = (POLYGON *) palloc(size);
SET_VARSIZE(poly, size);
poly->npts = 4;
poly->p[0].x = box->low.x;
poly->p[0].y = box->low.y;
poly->p[1].x = box->low.x;
poly->p[1].y = box->high.y;
poly->p[2].x = box->high.x;
poly->p[2].y = box->high.y;
poly->p[3].x = box->high.x;
poly->p[3].y = box->low.y;
box_fill(&poly->boundbox, box->high.x, box->low.x,
box->high.y, box->low.y);
PG_RETURN_POLYGON_P(poly);
}
Datum
poly_path(PG_FUNCTION_ARGS)
{
POLYGON *poly = PG_GETARG_POLYGON_P(0);
PATH *path;
int size;
int i;
size = offsetof(PATH, p[0]) +sizeof(path->p[0]) * poly->npts;
path = (PATH *) palloc(size);
SET_VARSIZE(path, size);
path->npts = poly->npts;
path->closed = TRUE;
for (i = 0; i < poly->npts; i++)
{
path->p[i].x = poly->p[i].x;
path->p[i].y = poly->p[i].y;
}
PG_RETURN_PATH_P(path);
}
/***********************************************************************
**
** Routines for circles.
**
***********************************************************************/
/*----------------------------------------------------------
* Formatting and conversion routines.
*---------------------------------------------------------*/
/* circle_in - convert a string to internal form.
*
* External format: (center and radius of circle)
* "((f8,f8)<f8>)"
* also supports quick entry style "(f8,f8,f8)"
*/
Datum
circle_in(PG_FUNCTION_ARGS)
{
char *str = PG_GETARG_CSTRING(0);
CIRCLE *circle;
char *s,
*cp;
int depth = 0;
circle = (CIRCLE *) palloc(sizeof(CIRCLE));
s = str;
while (isspace((unsigned char) *s))
s++;
if ((*s == LDELIM_C) || (*s == LDELIM))
{
depth++;
cp = (s + 1);
while (isspace((unsigned char) *cp))
cp++;
if (*cp == LDELIM)
s = cp;
}
if (!pair_decode(s, &circle->center.x, &circle->center.y, &s))
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type circle: \"%s\"", str)));
if (*s == DELIM)
s++;
while (isspace((unsigned char) *s))
s++;
if ((!single_decode(s, &circle->radius, &s)) || (circle->radius < 0))
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type circle: \"%s\"", str)));
while (depth > 0)
{
if ((*s == RDELIM)
|| ((*s == RDELIM_C) && (depth == 1)))
{
depth--;
s++;
while (isspace((unsigned char) *s))
s++;
}
else
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type circle: \"%s\"", str)));
}
if (*s != '\0')
ereport(ERROR,
(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
errmsg("invalid input syntax for type circle: \"%s\"", str)));
PG_RETURN_CIRCLE_P(circle);
}
/* circle_out - convert a circle to external form.
*/
Datum
circle_out(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
char *result;
char *cp;
result = palloc(2 * P_MAXLEN + 6);
cp = result;
*cp++ = LDELIM_C;
*cp++ = LDELIM;
if (!pair_encode(circle->center.x, circle->center.y, cp))
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("could not format \"circle\" value")));
cp += strlen(cp);
*cp++ = RDELIM;
*cp++ = DELIM;
if (!single_encode(circle->radius, cp))
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("could not format \"circle\" value")));
cp += strlen(cp);
*cp++ = RDELIM_C;
*cp = '\0';
PG_RETURN_CSTRING(result);
}
/*
* circle_recv - converts external binary format to circle
*/
Datum
circle_recv(PG_FUNCTION_ARGS)
{
StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
CIRCLE *circle;
circle = (CIRCLE *) palloc(sizeof(CIRCLE));
circle->center.x = pq_getmsgfloat8(buf);
circle->center.y = pq_getmsgfloat8(buf);
circle->radius = pq_getmsgfloat8(buf);
if (circle->radius < 0)
ereport(ERROR,
(errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
errmsg("invalid radius in external \"circle\" value")));
PG_RETURN_CIRCLE_P(circle);
}
/*
* circle_send - converts circle to binary format
*/
Datum
circle_send(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
StringInfoData buf;
pq_begintypsend(&buf);
pq_sendfloat8(&buf, circle->center.x);
pq_sendfloat8(&buf, circle->center.y);
pq_sendfloat8(&buf, circle->radius);
PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
}
/*----------------------------------------------------------
* Relational operators for CIRCLEs.
* <, >, <=, >=, and == are based on circle area.
*---------------------------------------------------------*/
/* circles identical?
*/
Datum
circle_same(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPeq(circle1->radius, circle2->radius) &&
FPeq(circle1->center.x, circle2->center.x) &&
FPeq(circle1->center.y, circle2->center.y));
}
/* circle_overlap - does circle1 overlap circle2?
*/
Datum
circle_overlap(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
circle1->radius + circle2->radius));
}
/* circle_overleft - is the right edge of circle1 at or left of
* the right edge of circle2?
*/
Datum
circle_overleft(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPle((circle1->center.x + circle1->radius),
(circle2->center.x + circle2->radius)));
}
/* circle_left - is circle1 strictly left of circle2?
*/
Datum
circle_left(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPlt((circle1->center.x + circle1->radius),
(circle2->center.x - circle2->radius)));
}
/* circle_right - is circle1 strictly right of circle2?
*/
Datum
circle_right(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPgt((circle1->center.x - circle1->radius),
(circle2->center.x + circle2->radius)));
}
/* circle_overright - is the left edge of circle1 at or right of
* the left edge of circle2?
*/
Datum
circle_overright(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPge((circle1->center.x - circle1->radius),
(circle2->center.x - circle2->radius)));
}
/* circle_contained - is circle1 contained by circle2?
*/
Datum
circle_contained(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPle((point_dt(&circle1->center, &circle2->center) + circle1->radius), circle2->radius));
}
/* circle_contain - does circle1 contain circle2?
*/
Datum
circle_contain(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPle((point_dt(&circle1->center, &circle2->center) + circle2->radius), circle1->radius));
}
/* circle_below - is circle1 strictly below circle2?
*/
Datum
circle_below(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPlt((circle1->center.y + circle1->radius),
(circle2->center.y - circle2->radius)));
}
/* circle_above - is circle1 strictly above circle2?
*/
Datum
circle_above(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPgt((circle1->center.y - circle1->radius),
(circle2->center.y + circle2->radius)));
}
/* circle_overbelow - is the upper edge of circle1 at or below
* the upper edge of circle2?
*/
Datum
circle_overbelow(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPle((circle1->center.y + circle1->radius),
(circle2->center.y + circle2->radius)));
}
/* circle_overabove - is the lower edge of circle1 at or above
* the lower edge of circle2?
*/
Datum
circle_overabove(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPge((circle1->center.y - circle1->radius),
(circle2->center.y - circle2->radius)));
}
/* circle_relop - is area(circle1) relop area(circle2), within
* our accuracy constraint?
*/
Datum
circle_eq(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPeq(circle_ar(circle1), circle_ar(circle2)));
}
Datum
circle_ne(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPne(circle_ar(circle1), circle_ar(circle2)));
}
Datum
circle_lt(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPlt(circle_ar(circle1), circle_ar(circle2)));
}
Datum
circle_gt(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPgt(circle_ar(circle1), circle_ar(circle2)));
}
Datum
circle_le(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPle(circle_ar(circle1), circle_ar(circle2)));
}
Datum
circle_ge(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
PG_RETURN_BOOL(FPge(circle_ar(circle1), circle_ar(circle2)));
}
/*----------------------------------------------------------
* "Arithmetic" operators on circles.
*---------------------------------------------------------*/
static CIRCLE *
circle_copy(CIRCLE *circle)
{
CIRCLE *result;
if (!PointerIsValid(circle))
return NULL;
result = (CIRCLE *) palloc(sizeof(CIRCLE));
memcpy((char *) result, (char *) circle, sizeof(CIRCLE));
return result;
}
/* circle_add_pt()
* Translation operator.
*/
Datum
circle_add_pt(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
Point *point = PG_GETARG_POINT_P(1);
CIRCLE *result;
result = circle_copy(circle);
result->center.x += point->x;
result->center.y += point->y;
PG_RETURN_CIRCLE_P(result);
}
Datum
circle_sub_pt(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
Point *point = PG_GETARG_POINT_P(1);
CIRCLE *result;
result = circle_copy(circle);
result->center.x -= point->x;
result->center.y -= point->y;
PG_RETURN_CIRCLE_P(result);
}
/* circle_mul_pt()
* Rotation and scaling operators.
*/
Datum
circle_mul_pt(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
Point *point = PG_GETARG_POINT_P(1);
CIRCLE *result;
Point *p;
result = circle_copy(circle);
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);
PG_RETURN_CIRCLE_P(result);
}
Datum
circle_div_pt(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
Point *point = PG_GETARG_POINT_P(1);
CIRCLE *result;
Point *p;
result = circle_copy(circle);
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);
PG_RETURN_CIRCLE_P(result);
}
/* circle_area - returns the area of the circle.
*/
Datum
circle_area(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
PG_RETURN_FLOAT8(circle_ar(circle));
}
/* circle_diameter - returns the diameter of the circle.
*/
Datum
circle_diameter(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
PG_RETURN_FLOAT8(2 * circle->radius);
}
/* circle_radius - returns the radius of the circle.
*/
Datum
circle_radius(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
PG_RETURN_FLOAT8(circle->radius);
}
/* circle_distance - returns the distance between
* two circles.
*/
Datum
circle_distance(PG_FUNCTION_ARGS)
{
CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
float8 result;
result = point_dt(&circle1->center, &circle2->center)
- (circle1->radius + circle2->radius);
if (result < 0)
result = 0;
PG_RETURN_FLOAT8(result);
}
Datum
circle_contain_pt(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
Point *point = PG_GETARG_POINT_P(1);
double d;
d = point_dt(&circle->center, point);
PG_RETURN_BOOL(d <= circle->radius);
}
Datum
pt_contained_circle(PG_FUNCTION_ARGS)
{
Point *point = PG_GETARG_POINT_P(0);
CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
double d;
d = point_dt(&circle->center, point);
PG_RETURN_BOOL(d <= circle->radius);
}
/* dist_pc - returns the distance between
* a point and a circle.
*/
Datum
dist_pc(PG_FUNCTION_ARGS)
{
Point *point = PG_GETARG_POINT_P(0);
CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
float8 result;
result = point_dt(point, &circle->center) - circle->radius;
if (result < 0)
result = 0;
PG_RETURN_FLOAT8(result);
}
/* circle_center - returns the center point of the circle.
*/
Datum
circle_center(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
Point *result;
result = (Point *) palloc(sizeof(Point));
result->x = circle->center.x;
result->y = circle->center.y;
PG_RETURN_POINT_P(result);
}
/* circle_ar - returns the area of the circle.
*/
static double
circle_ar(CIRCLE *circle)
{
return M_PI * (circle->radius * circle->radius);
}
/*----------------------------------------------------------
* Conversion operators.
*---------------------------------------------------------*/
Datum
cr_circle(PG_FUNCTION_ARGS)
{
Point *center = PG_GETARG_POINT_P(0);
float8 radius = PG_GETARG_FLOAT8(1);
CIRCLE *result;
result = (CIRCLE *) palloc(sizeof(CIRCLE));
result->center.x = center->x;
result->center.y = center->y;
result->radius = radius;
PG_RETURN_CIRCLE_P(result);
}
Datum
circle_box(PG_FUNCTION_ARGS)
{
CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
BOX *box;
double delta;
box = (BOX *) palloc(sizeof(BOX));
delta = circle->radius / sqrt(2.0);
box->high.x = circle->center.x + delta;
box->low.x = circle->center.x - delta;
box->high.y = circle->center.y + delta;
box->low.y = circle->center.y - delta;
PG_RETURN_BOX_P(box);
}
/* box_circle()
* Convert a box to a circle.
*/
Datum
box_circle(PG_FUNCTION_ARGS)
{
BOX *box = PG_GETARG_BOX_P(0);
CIRCLE *circle;
circle = (CIRCLE *) palloc(sizeof(CIRCLE));
circle->center.x = (box->high.x + box->low.x) / 2;
circle->center.y = (box->high.y + box->low.y) / 2;
circle->radius = point_dt(&circle->center, &box->high);
PG_RETURN_CIRCLE_P(circle);
}
Datum
circle_poly(PG_FUNCTION_ARGS)
{
int32 npts = PG_GETARG_INT32(0);
CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
POLYGON *poly;
int base_size,
size;
int i;
double angle;
double anglestep;
if (FPzero(circle->radius))
ereport(ERROR,
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
errmsg("cannot convert circle with radius zero to polygon")));
if (npts < 2)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("must request at least 2 points")));
base_size = sizeof(poly->p[0]) * npts;
size = offsetof(POLYGON, p[0]) +base_size;
/* Check for integer overflow */
if (base_size / npts != sizeof(poly->p[0]) || size <= base_size)
ereport(ERROR,
(errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
errmsg("too many points requested")));
poly = (POLYGON *) palloc0(size); /* zero any holes */
SET_VARSIZE(poly, size);
poly->npts = npts;
anglestep = (2.0 * M_PI) / npts;
for (i = 0; i < npts; i++)
{
angle = i * anglestep;
poly->p[i].x = circle->center.x - (circle->radius * cos(angle));
poly->p[i].y = circle->center.y + (circle->radius * sin(angle));
}
make_bound_box(poly);
PG_RETURN_POLYGON_P(poly);
}
/* poly_circle - convert polygon to circle
*
* 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)
{
POLYGON *poly = PG_GETARG_POLYGON_P(0);
CIRCLE *circle;
int i;
if (poly->npts < 2)
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("cannot convert empty polygon to circle")));
circle = (CIRCLE *) palloc(sizeof(CIRCLE));
circle->center.x = 0;
circle->center.y = 0;
circle->radius = 0;
for (i = 0; i < poly->npts; i++)
{
circle->center.x += poly->p[i].x;
circle->center.y += poly->p[i].y;
}
circle->center.x /= poly->npts;
circle->center.y /= poly->npts;
for (i = 0; i < poly->npts; i++)
circle->radius += point_dt(&poly->p[i], &circle->center);
circle->radius /= poly->npts;
if (FPzero(circle->radius))
ereport(ERROR,
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
errmsg("cannot convert empty polygon to circle")));
PG_RETURN_CIRCLE_P(circle);
}
/***********************************************************************
**
** Private routines for multiple types.
**
***********************************************************************/
/*
* Test to see if the point is inside the polygon, returns 1/0, or 2 if
* the point is on the polygon.
* Code adapted but not copied from integer-based routines in WN: A
* Server for the HTTP
* version 1.15.1, file wn/image.c
* http://hopf.math.northwestern.edu/index.html
* Description of algorithm: http://www.linuxjournal.com/article/2197
* http://www.linuxjournal.com/article/2029
*/
#define POINT_ON_POLYGON INT_MAX
static int
point_inside(Point *p, int npts, Point *plist)
{
double x0,
y0;
double prev_x,
prev_y;
int i = 0;
double x,
y;
int cross,
total_cross = 0;
if (npts <= 0)
return 0;
/* compute first polygon point relative to single point */
x0 = plist[0].x - p->x;
y0 = plist[0].y - p->y;
prev_x = x0;
prev_y = y0;
/* loop over polygon points and aggregate total_cross */
for (i = 1; i < npts; i++)
{
/* compute next polygon point relative to single point */
x = plist[i].x - p->x;
y = plist[i].y - p->y;
/* compute previous to current point crossing */
if ((cross = lseg_crossing(x, y, prev_x, prev_y)) == POINT_ON_POLYGON)
return 2;
total_cross += cross;
prev_x = x;
prev_y = y;
}
/* now do the first point */
if ((cross = lseg_crossing(x0, y0, prev_x, prev_y)) == POINT_ON_POLYGON)
return 2;
total_cross += cross;
if (total_cross != 0)
return 1;
return 0;
}
/* lseg_crossing()
* Returns +/-2 if line segment crosses the positive X-axis in a +/- direction.
* Returns +/-1 if one point is on the positive X-axis.
* Returns 0 if both points are on the positive X-axis, or there is no crossing.
* Returns POINT_ON_POLYGON if the segment contains (0,0).
* Wow, that is one confusing API, but it is used above, and when summed,
* can tell is if a point is in a polygon.
*/
static int
lseg_crossing(double x, double y, double prev_x, double prev_y)
{
double z;
int y_sign;
if (FPzero(y))
{ /* y == 0, on X axis */
if (FPzero(x)) /* (x,y) is (0,0)? */
return POINT_ON_POLYGON;
else if (FPgt(x, 0))
{ /* x > 0 */
if (FPzero(prev_y)) /* y and prev_y are zero */
/* prev_x > 0? */
return FPgt(prev_x, 0) ? 0 : POINT_ON_POLYGON;
return FPlt(prev_y, 0) ? 1 : -1;
}
else
{ /* x < 0, x not on positive X axis */
if (FPzero(prev_y))
/* prev_x < 0? */
return FPlt(prev_x, 0) ? 0 : POINT_ON_POLYGON;
return 0;
}
}
else
{ /* y != 0 */
/* compute y crossing direction from previous point */
y_sign = FPgt(y, 0) ? 1 : -1;
if (FPzero(prev_y))
/* previous point was on X axis, so new point is either off or on */
return FPlt(prev_x, 0) ? 0 : y_sign;
else if (FPgt(y_sign * prev_y, 0))
/* both above or below X axis */
return 0; /* same sign */
else
{ /* y and prev_y cross X-axis */
if (FPge(x, 0) && FPgt(prev_x, 0))
/* both non-negative so cross positive X-axis */
return 2 * y_sign;
if (FPlt(x, 0) && FPle(prev_x, 0))
/* both non-positive so do not cross positive X-axis */
return 0;
/* x and y cross axises, see URL above point_inside() */
z = (x - prev_x) * y - (y - prev_y) * x;
if (FPzero(z))
return POINT_ON_POLYGON;
return FPgt((y_sign * z), 0) ? 0 : 2 * y_sign;
}
}
}
static bool
plist_same(int npts, Point *p1, Point *p2)
{
int i,
ii,
j;
/* find match for first point */
for (i = 0; i < npts; i++)
{
if ((FPeq(p2[i].x, p1[0].x))
&& (FPeq(p2[i].y, p1[0].y)))
{
/* match found? then look forward through remaining points */
for (ii = 1, j = i + 1; ii < npts; ii++, j++)
{
if (j >= npts)
j = 0;
if ((!FPeq(p2[j].x, p1[ii].x))
|| (!FPeq(p2[j].y, p1[ii].y)))
{
#ifdef GEODEBUG
printf("plist_same- %d failed forward match with %d\n", j, ii);
#endif
break;
}
}
#ifdef GEODEBUG
printf("plist_same- ii = %d/%d after forward match\n", ii, npts);
#endif
if (ii == npts)
return TRUE;
/* match not found forwards? then look backwards */
for (ii = 1, j = i - 1; ii < npts; ii++, j--)
{
if (j < 0)
j = (npts - 1);
if ((!FPeq(p2[j].x, p1[ii].x))
|| (!FPeq(p2[j].y, p1[ii].y)))
{
#ifdef GEODEBUG
printf("plist_same- %d failed reverse match with %d\n", j, ii);
#endif
break;
}
}
#ifdef GEODEBUG
printf("plist_same- ii = %d/%d after reverse match\n", ii, npts);
#endif
if (ii == npts)
return TRUE;
}
}
return FALSE;
}
Reference in New Issue View Git Blame Copy Permalink