geo_ops.c

PostgreSQL 8.1.4的源码 适用于Linux下的开源数据库系统

 *				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)
	{							/* slope is infinite */
		m = (double) DBL_MAX;
	}
	else
		m = ((lseg->p[0].y - lseg->p[1].y) / (lseg->p[1].x - lseg->p[0].x));
	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 endpoints of the
	 * segment.
	 */

	/* intersection is on the line segment? */
	if ((ip = interpt_sl(lseg, ln)) != NULL)
	{
		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
	{
		/* intersection is not on line segment */
		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 [e-mail]; 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 ((d = 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 ((d = 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 = d = 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 th