geo_ops.c

关系型数据库 Postgresql 6.5.2

int4
path_npoints(PATH *path)
{
	if (!PointerIsValid(path))
		return 0;

	return path->npts;
}	/* path_npoints() */

PATH *
path_close(PATH *path)
{
	PATH	   *result;

	if (!PointerIsValid(path))
		return NULL;

	result = path_copy(path);
	result->closed = TRUE;

	return result;
}	/* path_close() */


PATH *
path_open(PATH *path)
{
	PATH	   *result;

	if (!PointerIsValid(path))
		return NULL;

	result = path_copy(path);
	result->closed = FALSE;

	return result;
}	/* path_open() */


static PATH *
path_copy(PATH *path)
{
	PATH	   *result;
	int			size;

	size = offsetof(PATH, p[0]) +(sizeof(path->p[0]) * path->npts);
	result = palloc(size);

	memmove((char *) result, (char *) path, size);
	return result;
}	/* path_copy() */


/* 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.
 */
bool
path_inter(PATH *p1, PATH *p2)
{
	BOX			b1,
				b2;
	int			i,
				j;
	LSEG		seg1,
				seg2;

	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_overlap(&b1, &b2))
		return FALSE;

	/* pairwise check lseg intersections */
	for (i = 0; i < p1->npts - 1; i++)
	{
		for (j = 0; j < p2->npts - 1; j++)
		{
			statlseg_construct(&seg1, &p1->p[i], &p1->p[i + 1]);
			statlseg_construct(&seg2, &p2->p[j], &p2->p[j + 1]);
			if (lseg_intersect(&seg1, &seg2))
				return TRUE;
		}
	}

	/* if we dropped through, no two segs intersected */
	return FALSE;
}	/* path_inter() */

/* path_distance()
 * This essentially does a cartesian product of the lsegs in the
 *	two paths, and finds the min distance between any two lsegs
 */
double *
path_distance(PATH *p1, PATH *p2)
{
	double	   *min = NULL,
			   *tmp;
	int			i,
				j;
	LSEG		seg1,
				seg2;

/*
	statlseg_construct(&seg1, &p1->p[0], &p1->p[1]);
	statlseg_construct(&seg2, &p2->p[0], &p2->p[1]);
	min = lseg_distance(&seg1, &seg2);
*/

	for (i = 0; i < p1->npts - 1; i++)
		for (j = 0; j < p2->npts - 1; j++)
		{
			statlseg_construct(&seg1, &p1->p[i], &p1->p[i + 1]);
			statlseg_construct(&seg2, &p2->p[j], &p2->p[j + 1]);

			tmp = lseg_distance(&seg1, &seg2);
			if ((min == NULL) || (*min < *tmp))
			{
				if (min != NULL)
					pfree(min);
				min = tmp;
			}
			else
				pfree(tmp);
		}

	return min;
}	/* path_distance() */


/*----------------------------------------------------------
 *	"Arithmetic" operations.
 *---------------------------------------------------------*/

double *
path_length(PATH *path)
{
	double	   *result;
	int			i;

	result = palloc(sizeof(double));

	*result = 0;
	for (i = 0; i < (path->npts - 1); i++)
		*result += point_dt(&path->p[i], &path->p[i + 1]);

	return result;
}	/* path_length() */


#ifdef NOT_USED
double
path_ln(PATH *path)
{
	double		result;
	int			i;

	result = 0;
	for (i = 0; i < (path->npts - 1); i++)
		result += point_dt(&path->p[i], &path->p[i + 1]);

	return result;
}	/* path_ln() */

#endif

/***********************************************************************
 **
 **		Routines for 2D points.
 **
 ***********************************************************************/

/*----------------------------------------------------------
 *	String to point, point to string conversion.
 *		External format:
 *				"(x,y)"
 *				"x,y"
 *---------------------------------------------------------*/

Point *
point_in(char *str)
{
	Point	   *point;

	double		x,
				y;
	char	   *s;

	if (!PointerIsValid(str))
		elog(ERROR, "Bad (null) point external representation");

	if (!pair_decode(str, &x, &y, &s) || (strlen(s) > 0))
		elog(ERROR, "Bad point external representation '%s'", str);

	point = palloc(sizeof(Point));

	point->x = x;
	point->y = y;

	return point;
}	/* point_in() */

char *
point_out(Point *pt)
{
	if (!PointerIsValid(pt))
		return NULL;

	return path_encode(-1, 1, pt);
}	/* point_out() */


static Point *
point_construct(double x, double y)
{
	Point	   *result = 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 = 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).
 *---------------------------------------------------------*/

bool
point_left(Point *pt1, Point *pt2)
{
	return FPlt(pt1->x, pt2->x);
}

bool
point_right(Point *pt1, Point *pt2)
{
	return FPgt(pt1->x, pt2->x);
}

bool
point_above(Point *pt1, Point *pt2)
{
	return FPgt(pt1->y, pt2->y);
}

bool
point_below(Point *pt1, Point *pt2)
{
	return FPlt(pt1->y, pt2->y);
}

bool
point_vert(Point *pt1, Point *pt2)
{
	return FPeq(pt1->x, pt2->x);
}

bool
point_horiz(Point *pt1, Point *pt2)
{
	return FPeq(pt1->y, pt2->y);
}

bool
point_eq(Point *pt1, Point *pt2)
{
	return point_horiz(pt1, pt2) && point_vert(pt1, pt2);
}

bool
point_ne(Point *pt1, Point *pt2)
{
	return !point_eq(pt1, pt2);
}

/*----------------------------------------------------------
 *	"Arithmetic" operators on points.
 *---------------------------------------------------------*/

int32
pointdist(Point *p1, Point *p2)
{
	int32		result;

	result = point_dt(p1, p2);
	return result;
}

double *
point_distance(Point *pt1, Point *pt2)
{
	double	   *result = palloc(sizeof(double));

	*result = HYPOT(pt1->x - pt2->x, pt1->y - pt2->y);
	return result;
}


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);
}

double *
point_slope(Point *pt1, Point *pt2)
{
	double	   *result = palloc(sizeof(double));

	if (point_vert(pt1, pt2))
		*result = (double) DBL_MAX;
	else
		*result = (pt1->y - pt2->y) / (pt1->x - pt1->x);
	return result;
}


double
point_sl(Point *pt1, Point *pt2)
{
	return (point_vert(pt1, pt2)
			? (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)"
 *---------------------------------------------------------*/

LSEG *
lseg_in(char *str)
{
	LSEG	   *lseg;

	int			isopen;
	char	   *s;

	if (!PointerIsValid(str))
		elog(ERROR, " Bad (null) lseg external representation", NULL);

	lseg = palloc(sizeof(LSEG));

	if ((!path_decode(TRUE, 2, str, &isopen, &s, &(lseg->p[0])))
		|| (*s != '\0'))
		elog(ERROR, "Bad lseg external representation '%s'", str);

#ifdef NOT_USED
	lseg->m = point_sl(&lseg->p[0], &lseg->p[1]);
#endif

	return lseg;
}	/* lseg_in() */


char *
lseg_out(LSEG *ls)
{
	if (!PointerIsValid(ls))
		return NULL;

	return path_encode(FALSE, 2, (Point *) &(ls->p[0]));
}	/* lseg_out() */


/* lseg_construct -
 *		form a LSEG from two Points.
 */
LSEG *
lseg_construct(Point *pt1, Point *pt2)
{
	LSEG	   *result = 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

	return 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
}

double *
lseg_length(LSEG *lseg)
{
	double	   *result;

	if (!PointerIsValid(lseg))
		return NULL;

	result = point_distance(&lseg->p[0], &lseg->p[1]);

	return result;
}	/* lseg_length() */

/*----------------------------------------------------------
 *	Relative position routines.
 *---------------------------------------------------------*/

/*
 **  find intersection of the two lines, and see if it falls on
 **  both segments.
 */
bool
lseg_intersect(LSEG *l1, LSEG *l2)
{
	LINE	   *ln;
	Point	   *interpt;
	bool		retval;

	ln = line_construct_pp(&l2->p[0], &l2->p[1]);
	interpt = interpt_sl(l1, ln);

	if (interpt != NULL && on_ps(interpt, l2))	/* interpt on l1 and l2 */
		retval = TRUE;
	else
		retval = FALSE;
	if (interpt != NULL)
		pfree(interpt);
	pfree(ln);
	return retval;
}

bool
lseg_parallel(LSEG *l1, LSEG *l2)
{
#ifdef NOT_USED
	return FPeq(l1->m, l2->m);
#endif
	return (FPeq(point_sl(&(l1->p[0]), &(l1->p[1])),
				 point_sl(&(l2->p[0]), &(l2->p[1]))));
}	/* lseg_parallel() */

/* 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
 */
bool
lseg_perp(LSEG *l1, LSEG *l2)
{
	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))
		return FPeq(m2, DBL_MAX);
	else if (FPzero(m2))
		return FPeq(m1, DBL_MAX);

	return FPeq(m1 / m2, -1.0);
}	/* lseg_perp() */

bool
lseg_vertical(LSEG *lseg)
{
	return FPeq(lseg->p[0].x, lseg->p[1].x);
}

bool
lseg_horizontal(LSEG *lseg)
{
	return FPeq(lseg->p[0].y, lseg->p[1].y);
}


bool
lseg_eq(LSEG *l1, LSEG *l2)
{
	return (FPeq(l1->p[0].x, l2->p[0].x) &&
			FPeq(l1->p[1].y, l2->p[1].y) &&
			FPeq(l1->p[0].x, l2->p[0].x) &&
			FPeq(l1->p[1].y, l2->p[1].y));
}	/* lseg_eq() */

bool
lseg_ne(LSEG *l1, LSEG *l2)
{
	return (!FPeq(l1->p[0].x, l2->p[0].x) ||
			!FPeq(l1->p[1].y, l2->p[1].y) ||
			!FPeq(l1->p[0].x, l2->p[0].x) ||
			!FPeq(l1->p[1].y, l2->p[1].y));
}	/* lseg_ne() */

bool
lseg_lt(LSEG *l1, LSEG *l2)
{
	return FPlt(point_dt(&l1->p[0], &l1->p[1]), point_dt(&l2->p[0], &l2->p[1]));
}	/* lseg_lt() */

bool
lseg_le(LSEG *l1, LSEG *l2)
{
	return FPle(point_dt(&l1->p[0], &l1->p[1]), point_dt(&l2->p[0], &l2->p[1]));
}	/* lseg_le() */

bool
lseg_gt(LSEG *l1, LSEG *l2)
{
	return FPgt(point_dt(&l1->p[0], &l1->p[1]), point_dt(&l2->p[0], &l2->p[1]));
}	/* lseg_gt() */

bool
lseg_ge(LSEG *l1, LSEG *l2)
{
	return FPge(point_dt(&l1->p[0], &l1->p[1]), point_dt(&l2->p[0], &l2->p[1]));
}	/* lseg_ge() */


/*----------------------------------------------------------
 *	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.
 */
double *
lseg_distance(LSEG *l1, LSEG *l2)
{
	double	   *result = palloc(sizeof(double));

	*result = lseg_dt(l1, l2);

	return result;
}

/* 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	   *d,
				result;

	if (lseg_intersect(l1, l2))
		return 0.0;

	d = dist_ps(&l1->p[0], l2);
	result = *d;
	pfree(d);
	d = dist_ps(&l1->p[1], l2);
	result = Min(result, *d);