geo_ops.c

关系型数据库 Postgresql 6.5.2

/*		box_distance	-		returns the distance between the
 *								  center points of two boxes.
 */
double *
box_distance(BOX *box1, BOX *box2)
{
	double	   *result = palloc(sizeof(double));
	Point	   *a,
			   *b;

	a = box_center(box1);
	b = box_center(box2);
	*result = HYPOT(a->x - b->x, a->y - b->y);

	pfree(a);
	pfree(b);
	return result;
}


/*		box_center		-		returns the center point of the box.
 */
Point *
box_center(BOX *box)
{
	Point	   *result = palloc(sizeof(Point));

	result->x = (box->high.x + box->low.x) / 2.0;
	result->y = (box->high.y + box->low.y) / 2.0;

	return result;
}


/*		box_ar	-		returns the area of the box.
 */
static double
box_ar(BOX *box)
{
	return box_wd(box) * box_ht(box);
}


/*		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;
}


/*		box_dt	-		returns the distance between the
 *						  center points of two boxes.
 */
#ifdef NOT_USED
static double
box_dt(BOX *box1, BOX *box2)
{
	double		result;
	Point	   *a,
			   *b;

	a = box_center(box1);
	b = box_center(box2);
	result = HYPOT(a->x - b->x, a->y - b->y);

	pfree(a);
	pfree(b);
	return result;
}

#endif

/*----------------------------------------------------------
 *	Funky operations.
 *---------------------------------------------------------*/

/*		box_intersect	-
 *				returns the overlapping portion of two boxes,
 *				  or NULL if they do not intersect.
 */
BOX *
box_intersect(BOX *box1, BOX *box2)
{
	BOX		   *result;

	if (!box_overlap(box1, box2))
		return NULL;

	result = 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);

	return result;
}


/*		box_diagonal	-
 *				returns a line segment which happens to be the
 *				  positive-slope diagonal of "box".
 *				provided, of course, we have LSEGs.
 */
LSEG *
box_diagonal(BOX *box)
{
	Point		p1,
				p2;

	p1.x = box->high.x;
	p1.y = box->high.y;
	p2.x = box->low.x;
	p2.y = box->low.y;
	return lseg_construct(&p1, &p2);
}

/***********************************************************************
 **
 **		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.
 **
 ***********************************************************************/

LINE *
line_in(char *str)
{
	LINE	   *line;

#ifdef ENABLE_LINE_TYPE
	LSEG		lseg;
	int			isopen;
	char	   *s;

#endif

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

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

	line = line_construct_pp(&(lseg.p[0]), &(lseg.p[1]));
#else
	elog(ERROR, "line not yet implemented");
	line = NULL;
#endif

	return line;
}	/* line_in() */


char *
line_out(LINE *line)
{
	char	   *result;

#ifdef ENABLE_LINE_TYPE
	LSEG		lseg;

#endif

	if (!PointerIsValid(line))
		return NULL;

#ifdef ENABLE_LINE_TYPE
	if (FPzero(line->B))
	{							/* vertical */
		/* use "x = C" */
		result->A = -1;
		result->B = 0;
		result->C = pt1->x;
#ifdef GEODEBUG
		printf("line_construct_pp- 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_construct_pp- line is horizontal\n");
#endif
#ifdef NOT_USED
		result->m = 0.0;
#endif

	}
	else
	{
	}

	if (line_horizontal(line))
	{
	}
	else if (line_vertical(line))
	{
	}
	else
	{
	}

	return path_encode(TRUE, 2, (Point *) &(ls->p[0]));
#else
	elog(ERROR, "line not yet implemented");
	result = NULL;
#endif

	return result;
}	/* line_out() */


/*----------------------------------------------------------
 *	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 = palloc(sizeof(LINE));

	/* use "mx - y + yinter = 0" */
	result->A = m;
	result->B = -1.0;
	result->C = pt->y - m * pt->x;

#ifdef NOT_USED
	result->m = m;
#endif

	return result;
}	/* line_construct_pm() */


/* line_construct_pp()
 * two points
 */
LINE *
line_construct_pp(Point *pt1, Point *pt2)
{
	LINE	   *result = palloc(sizeof(LINE));

	if (FPeq(pt1->x, pt2->x))
	{							/* vertical */
		/* use "x = C" */
		result->A = -1;
		result->B = 0;
		result->C = pt1->x;
#ifdef GEODEBUG
		printf("line_construct_pp- line is vertical\n");
#endif
#ifdef NOT_USED
		result->m = DBL_MAX;
#endif

	}
	else if (FPeq(pt1->y, pt2->y))
	{							/* horizontal */
		/* use "x = C" */
		result->A = 0;
		result->B = -1;
		result->C = pt1->y;
#ifdef GEODEBUG
		printf("line_construct_pp- line is horizontal\n");
#endif
#ifdef NOT_USED
		result->m = 0.0;
#endif

	}
	else
	{
		/* use "mx - y + yinter = 0" */
#ifdef NOT_USED
		result->A = (pt1->y - pt2->y) / (pt1->x - pt2->x);
#endif
		result->A = (pt2->y - pt1->y) / (pt2->x - pt1->x);
		result->B = -1.0;
		result->C = pt1->y - result->A * pt1->x;
#ifdef GEODEBUG
		printf("line_construct_pp- line is neither vertical nor horizontal (diffs x=%.*g, y=%.*g\n",
			   digits8, (pt2->x - pt1->x), digits8, (pt2->y - pt1->y));
#endif
#ifdef NOT_USED
		result->m = result->A;
#endif
	}
	return result;
}	/* line_construct_pp() */


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

bool
line_intersect(LINE *l1, LINE *l2)
{
	return !line_parallel(l1, l2);
}

bool
line_parallel(LINE *l1, LINE *l2)
{
#ifdef NOT_USED
	return FPeq(l1->m, l2->m);
#endif
	if (FPzero(l1->B))
		return FPzero(l2->B);

	return FPeq(l2->A, l1->A * (l2->B / l1->B));
}	/* line_parallel() */

bool
line_perp(LINE *l1, LINE *l2)
{
#ifdef NOT_USED
	if (l1->m)
		return FPeq(l2->m / l1->m, -1.0);
	else if (l2->m)
		return FPeq(l1->m / l2->m, -1.0);
#endif
	if (FPzero(l1->A))
		return FPzero(l2->B);
	else if (FPzero(l1->B))
		return FPzero(l2->A);

	return FPeq(((l1->A * l2->B) / (l1->B * l2->A)), -1.0);
}	/* line_perp() */

bool
line_vertical(LINE *line)
{
#ifdef NOT_USED
	return FPeq(line->A, -1.0) && FPzero(line->B);
#endif
	return FPzero(line->B);
}	/* line_vertical() */

bool
line_horizontal(LINE *line)
{
#ifdef NOT_USED
	return FPzero(line->m);
#endif
	return FPzero(line->A);
}	/* line_horizontal() */

bool
line_eq(LINE *l1, LINE *l2)
{
	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;

	return (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.
 */
double *
line_distance(LINE *l1, LINE *l2)
{
	double	   *result = palloc(sizeof(double));
	Point	   *tmp;

	if (line_intersect(l1, l2))
	{
		*result = 0.0;
		return result;
	}
	if (line_vertical(l1))
		*result = fabs(l1->C - l2->C);
	else
	{
		tmp = point_construct(0.0, l1->C);
		result = dist_pl(tmp, l2);
		pfree(tmp);
	}
	return result;
}

/* line_interpt()
 * Point where two lines l1, l2 intersect (if any)
 */
Point *
line_interpt(LINE *l1, LINE *l2)
{
	Point	   *result;
	double		x,
				y;

	if (line_parallel(l1, l2))
		return NULL;
#ifdef NOT_USED
	if (line_vertical(l1))
		result = point_construct(l2->m * l1->C + l2->C, l1->C);
	else if (line_vertical(l2))
		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 (line_vertical(l1))
	{
#ifdef NOT_USED
		x = l1->C;
		y = -((l2->A * x + l2->C) / l2->B);
#endif
		x = l1->C;
		y = (l2->A * x + l2->C);

	}
	else if (line_vertical(l2))
	{
#ifdef NOT_USED
		x = l2->C;
		y = -((l1->A * x + l1->C) / l1->B);
#endif
		x = l2->C;
		y = (l1->A * x + l1->C);

	}
	else
	{
#ifdef NOT_USED
		x = (l2->B * l1->C - l1->B * l2->C) / (l2->A * l1->B - l1->A * l2->B);
		y = -((l1->A * x + l1->C) / l1->B);
#endif
		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",
		   digits8, l1->A, digits8, l1->B, digits8, l1->C, digits8, l2->A, digits8, l2->B, digits8, l2->C);
	printf("line_interpt- lines intersect at (%.*g,%.*g)\n", digits8, x, digits8, y);
#endif
	return result;
}	/* line_interpt() */


/***********************************************************************
 **
 **		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,... )"
 *---------------------------------------------------------*/

PATH *
path_in(char *str)
{
	PATH	   *path;

	int			isopen;
	char	   *s;
	int			npts;
	int			size;
	int			depth = 0;

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

	if ((npts = pair_count(str, ',')) <= 0)
		elog(ERROR, "Bad path external representation '%s'", str);

	s = str;
	while (isspace(*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 = palloc(size);

	path->size = size;
	path->npts = npts;

	if ((!path_decode(TRUE, npts, s, &isopen, &s, &(path->p[0])))
		&& (!((depth == 0) && (*s == '\0'))) && !((depth >= 1) && (*s == RDELIM)))
		elog(ERROR, "Bad path external representation '%s'", str);

	path->closed = (!isopen);

	return path;
}	/* path_in() */


char *
path_out(PATH *path)
{
	if (!PointerIsValid(path))
		return NULL;

	return path_encode(path->closed, path->npts, (Point *) &(path->p[0]));
}	/* path_out() */


/*----------------------------------------------------------
 *	Relational operators.
 *		These are based on the path cardinality,
 *		as stupid as that sounds.
 *
 *		Better relops and access methods coming soon.
 *---------------------------------------------------------*/

bool
path_n_lt(PATH *p1, PATH *p2)
{
	return (p1->npts < p2->npts);
}

bool
path_n_gt(PATH *p1, PATH *p2)
{
	return (p1->npts > p2->npts);
}

bool
path_n_eq(PATH *p1, PATH *p2)
{
	return (p1->npts == p2->npts);
}

bool
path_n_le(PATH *p1, PATH *p2)
{
	return (p1->npts <= p2->npts);
}

bool
path_n_ge(PATH *p1, PATH *p2)
{
	return (p1->npts >= p2->npts);
}


/*----------------------------------------------------------
 * Conversion operators.
 *---------------------------------------------------------*/

bool
path_isclosed(PATH *path)
{
	if (!PointerIsValid(path))
		return FALSE;

	return path->closed;
}	/* path_isclosed() */

bool
path_isopen(PATH *path)
{
	if (!PointerIsValid(path))
		return FALSE;

	return !path->closed;
}	/* path_isopen() */