geo_ops.c

关系型数据库 Postgresql 6.5.2

	LSEG		lseg,
				seg;
	Point		point;
	double		dist,
			   *d;

	if (on_pb(pt, box))
		return pt;

	/* pairwise check lseg distances */
	point.x = box->low.x;
	point.y = box->high.y;
	statlseg_construct(&lseg, &box->low, &point);
	dist = *(d = dist_ps(pt, &lseg));
	pfree(d);

	statlseg_construct(&seg, &box->high, &point);
	if (*(d = dist_ps(pt, &seg)) < dist)
	{
		dist = *d;
		memcpy(&lseg, &seg, sizeof(lseg));
	}
	pfree(d);

	point.x = box->high.x;
	point.y = box->low.y;
	statlseg_construct(&seg, &box->low, &point);
	if (*(d = dist_ps(pt, &seg)) < dist)
	{
		dist = *d;
		memcpy(&lseg, &seg, sizeof(lseg));
	}
	pfree(d);

	statlseg_construct(&seg, &box->high, &point);
	if (*(d = dist_ps(pt, &seg)) < dist)
	{
		dist = *d;
		memcpy(&lseg, &seg, sizeof(lseg));
	}
	pfree(d);

	return close_ps(pt, &lseg);
}	/* close_pb() */

/* 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
 */
Point *
close_sl(LSEG *lseg, LINE *line)
{
	Point	   *result;
	double	   *d1,
			   *d2;

	result = interpt_sl(lseg, line);
	if (result)
		return result;

	d1 = dist_pl(&lseg->p[0], line);
	d2 = dist_pl(&lseg->p[1], line);
	if (d1 < d2)
		result = point_copy(&lseg->p[0]);
	else
		result = point_copy(&lseg->p[1]);

	pfree(d1);
	pfree(d2);
	return result;
}

/* close_ls()
 * Closest point on line segment to line.
 */
Point *
close_ls(LINE *line, LSEG *lseg)
{
	Point	   *result;
	double	   *d1,
			   *d2;

	result = interpt_sl(lseg, line);
	if (result)
		return result;

	d1 = dist_pl(&lseg->p[0], line);
	d2 = dist_pl(&lseg->p[1], line);
	if (d1 < d2)
		result = point_copy(&lseg->p[0]);
	else
		result = point_copy(&lseg->p[1]);

	pfree(d1);
	pfree(d2);
	return result;
}	/* close_ls() */

/* close_sb()
 * Closest point on or in box to line segment.
 */
Point *
close_sb(LSEG *lseg, BOX *box)
{
	Point	   *result;
	Point	   *pt;
	Point		point;

	LSEG		bseg,
				seg;
	double		dist,
				d;

	/* segment intersects box? then just return closest point to center */
	if (inter_sb(lseg, box))
	{
		pt = box_center(box);
		result = close_ps(pt, lseg);
		pfree(pt);
		return result;
	}

	/* 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 */
	return close_lseg(lseg, &bseg);
}	/* close_sb() */

Point *
close_lb(LINE *line, BOX *box)
{
	/* think about this one for a while */
	elog(ERROR, "close_lb not implemented", NULL);

	return NULL;
}

/*---------------------------------------------------------------------
 *		on_
 *				Whether one object lies completely within another.
 *-------------------------------------------------------------------*/

/* on_pl -
 *		Does the point satisfy the equation?
 */
bool
on_pl(Point *pt, LINE *line)
{
	if (!PointerIsValid(pt) || !PointerIsValid(line))
		return FALSE;

	return 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
 */
bool
on_ps(Point *pt, LSEG *lseg)
{
	if (!PointerIsValid(pt) || !PointerIsValid(lseg))
		return FALSE;

	return (FPeq(point_dt(pt, &lseg->p[0]) + point_dt(pt, &lseg->p[1]),
				 point_dt(&lseg->p[0], &lseg->p[1])));
}

bool
on_pb(Point *pt, BOX *box)
{
	if (!PointerIsValid(pt) || !PointerIsValid(box))
		return FALSE;

	return (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)
 */
#define NEXT(A) ((A+1) % path->npts)	/* cyclic "i+1" */

bool
on_ppath(Point *pt, PATH *path)
{
	int			i,
				n;
	double		a,
				b;

	if (!PointerIsValid(pt) || !PointerIsValid(path))
		return FALSE;

	/*-- 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])))
				return TRUE;
			a = b;
		}
		return FALSE;
	}

	/*-- CLOSED --*/
	return point_inside(pt, path->npts, path->p);
}	/* on_ppath() */


bool
on_sl(LSEG *lseg, LINE *line)
{
	if (!PointerIsValid(lseg) || !PointerIsValid(line))
		return FALSE;

	return on_pl(&lseg->p[0], line) && on_pl(&lseg->p[1], line);
}	/* on_sl() */

bool
on_sb(LSEG *lseg, BOX *box)
{
	if (!PointerIsValid(lseg) || !PointerIsValid(box))
		return FALSE;

	return on_pb(&lseg->p[0], box) && on_pb(&lseg->p[1], box);
}	/* on_sb() */

/*---------------------------------------------------------------------
 *		inter_
 *				Whether one object intersects another.
 *-------------------------------------------------------------------*/

bool
inter_sl(LSEG *lseg, LINE *line)
{
	Point	   *tmp;

	if (!PointerIsValid(lseg) || !PointerIsValid(line))
		return FALSE;

	tmp = interpt_sl(lseg, line);
	if (tmp)
	{
		pfree(tmp);
		return TRUE;
	}
	return FALSE;
}

/* 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
 */
bool
inter_sb(LSEG *lseg, BOX *box)
{
	BOX			lbox;
	LSEG		bseg;
	Point		point;

	if (!PointerIsValid(lseg) || !PointerIsValid(box))
		return FALSE;

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

	/* an endpoint of segment is inside box? then clearly intersects */
	if (on_pb(&lseg->p[0], box) || on_pb(&lseg->p[1], box))
		return TRUE;

	/* pairwise check lseg intersections */
	point.x = box->low.x;
	point.y = box->high.y;
	statlseg_construct(&bseg, &box->low, &point);
	if (lseg_intersect(&bseg, lseg))
		return TRUE;

	statlseg_construct(&bseg, &box->high, &point);
	if (lseg_intersect(&bseg, lseg))
		return TRUE;

	point.x = box->high.x;
	point.y = box->low.y;
	statlseg_construct(&bseg, &box->low, &point);
	if (lseg_intersect(&bseg, lseg))
		return TRUE;

	statlseg_construct(&bseg, &box->high, &point);
	if (lseg_intersect(&bseg, lseg))
		return TRUE;

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

/* inter_lb()
 * Do line and box intersect?
 */
bool
inter_lb(LINE *line, BOX *box)
{
	LSEG		bseg;
	Point		p1,
				p2;

	if (!PointerIsValid(line) || !PointerIsValid(box))
		return FALSE;

	/* 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 (inter_sl(&bseg, line))
		return TRUE;
	p1.x = box->high.x;
	p1.y = box->high.y;
	statlseg_construct(&bseg, &p1, &p2);
	if (inter_sl(&bseg, line))
		return TRUE;
	p2.x = box->high.x;
	p2.y = box->low.y;
	statlseg_construct(&bseg, &p1, &p2);
	if (inter_sl(&bseg, line))
		return TRUE;
	p1.x = box->low.x;
	p1.y = box->low.y;
	statlseg_construct(&bseg, &p1, &p2);
	if (inter_sl(&bseg, line))
		return TRUE;

	/* if we dropped through, no intersection */
	return 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
		elog(ERROR, "Unable to create bounding box for empty polygon", NULL);
}

/*------------------------------------------------------------------
 * 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)"
 *------------------------------------------------------------------*/
POLYGON    *
poly_in(char *str)
{
	POLYGON    *poly;
	int			npts;
	int			size;
	int			isopen;
	char	   *s;

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

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

	size = offsetof(POLYGON, p[0]) +(sizeof(poly->p[0]) * npts);
	poly = palloc(size);

	MemSet((char *) poly, 0, size);		/* zero any holes */
	poly->size = size;
	poly->npts = npts;

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

	make_bound_box(poly);

	return poly;
}	/* poly_in() */

/*---------------------------------------------------------------
 * poly_out - convert internal POLYGON representation to the
 *			  character string format "((f8,f8),...,(f8,f8))"
 *			  also support old format "(f8,f8,...,f8,f8)"
 *---------------------------------------------------------------*/
char *
poly_out(POLYGON *poly)
{
	if (!PointerIsValid(poly))
		return NULL;

	return path_encode(TRUE, poly->npts, &(poly->p[0]));
}	/* poly_out() */


/*-------------------------------------------------------
 * 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?
 *-------------------------------------------------------*/
bool
poly_left(POLYGON *polya, POLYGON *polyb)
{
	return polya->boundbox.high.x < polyb->boundbox.low.x;
}

/*-------------------------------------------------------
 * Is polygon A overlapping or left of polygon B? i.e. is
 * the left most point of A left of the right most point
 * of B?
 *-------------------------------------------------------*/
bool
poly_overleft(POLYGON *polya, POLYGON *polyb)
{
	return polya->boundbox.low.x <= polyb->boundbox.high.x;
}

/*-------------------------------------------------------
 * 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?
 *-------------------------------------------------------*/
bool
poly_right(POLYGON *polya, POLYGON *polyb)
{
	return polya->boundbox.low.x > polyb->boundbox.high.x;
}

/*-------------------------------------------------------
 * Is polygon A overlapping or right of polygon B? i.e. is
 * the right most point of A right of the left most point
 * of B?
 *-------------------------------------------------------*/
bool
poly_overright(POLYGON *polya, POLYGON *polyb)
{
	return polya->boundbox.high.x > polyb->boundbox.low.x;
}

/*-------------------------------------------------------
 * 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.
 *-------------------------------------------------------*/
bool
poly_same(POLYGON *polya, POLYGON *polyb)
{
	if (!PointerIsValid(polya) || !PointerIsValid(polyb))
		return FALSE;

	if (polya->npts != polyb->npts)
		return FALSE;

	return plist_same(polya->npts, polya->p, polyb->p);

#ifdef NOT_USED
	for (i = 0; i < polya->npts; i++)
	{
		if ((polya->p[i].x != polyb->p[i].x)
			|| (polya->p[i].y != polyb->p[i].y))
			return FALSE;
	}
	return TRUE;
#endif
}	/* poly_same() */

/*-----------------------------------------------------------------
 * Determine if polygon A overlaps polygon B by determining if
 * their bounding boxes overlap.
 *-----------------------------------------------------------------*/
bool
poly_overlap(POLYGON *polya, POLYGON *polyb)
{
	return box_overlap(&(polya->boundbox), &(polyb->boundbox));
}


/*-----------------------------------------------------------------
 * Determine if polygon A contains polygon B by determining if A's
 * bounding box contains B's bounding box.
 *-----------------------------------------------------------------*/
#ifdef NOT_USED
bool
poly_contain(POLYGON *polya, POLYGON *polyb)
{
	return box_contain(&(polya->boundbox), &(polyb->boundbox));
}

#endif

bool
poly_contain(POLYGON *polya, POLYGON *polyb)
{
	int			i;

	if (!PointerIsValid(polya) || !PointerIsValid(polyb))
		return FALSE;

	if (box_contain(&(polya->boundbox), &(polyb->boundbox)))
	{
		for (i = 0; i < polyb->npts; i++)
		{
			if (point_inside(&(polyb->p[i]), polya->npts, &(polya->p[0])) == 0)
			{
#if GEODEBUG
				printf("poly_contain- point (%f,%f) not in polygon\n", polyb->p[i].x, polyb->p[i].y);
#endif
				return FALSE;
			}
		}
		for (i = 0; i < polya->npts; i++)
		{
			if (point_inside(&(polya->p[i]), polyb->npts, &(polyb->p[0])) == 1)
			{
#if GEODEBUG
				printf("poly_contain- point (%f,%f) in polygon\n", polya->p[i].x, polya->p[i].y);
#endif
				return FALSE;
			}
		}
		return TRUE;
	}
#if GEODEBUG
	printf("poly_contain- bound box ((%f,%f),(%f,%f)) not insid