geo_ops.c

关系型数据库 Postgresql 6.5.2

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

	return result;
}	/* lseg_dt() */


Point *
lseg_center(LSEG *lseg)
{
	Point	   *result;

	if (!PointerIsValid(lseg))
		return NULL;

	result = palloc(sizeof(Point));

	result->x = (lseg->p[0].x - lseg->p[1].x) / 2;
	result->y = (lseg->p[0].y - lseg->p[1].y) / 2;

	return result;
}	/* lseg_center() */


/* lseg_interpt -
 *		Find the intersection point of two segments (if any).
 *		Find the intersection of the appropriate lines; if the
 *		point is not on a given segment, there is no valid segment
 *		intersection point at all.
 * 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
 */
Point *
lseg_interpt(LSEG *l1, LSEG *l2)
{
	Point	   *result;
	LINE	   *tmp1,
			   *tmp2;

	if (!PointerIsValid(l1) || !PointerIsValid(l2))
		return NULL;

	tmp1 = line_construct_pp(&l1->p[0], &l1->p[1]);
	tmp2 = line_construct_pp(&l2->p[0], &l2->p[1]);
	result = line_interpt(tmp1, tmp2);
	if (PointerIsValid(result))
	{
		if (on_ps(result, l1))
		{
			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;
			}
		}
		else
		{
			pfree(result);
			result = NULL;
		}
	}
	pfree(tmp1);
	pfree(tmp2);

	return result;
}	/* lseg_interpt() */

/***********************************************************************
 **
 **		Routines for position comparisons of differently-typed
 **				2D objects.
 **
 ***********************************************************************/

#define ABOVE	1
#define BELOW	0
#define UNDEF	-1


/*---------------------------------------------------------------------
 *		dist_
 *				Minimum distance from one object to another.
 *-------------------------------------------------------------------*/

double *
dist_pl(Point *pt, LINE *line)
{
	double	   *result = palloc(sizeof(double));

	*result = (line->A * pt->x + line->B * pt->y + line->C) /
		HYPOT(line->A, line->B);

	return result;
}

double *
dist_ps(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
	{
#ifdef NOT_USED
		m = (-1) * (lseg->p[1].y - lseg->p[0].y) /
			(lseg->p[1].x - lseg->p[0].x);
#endif
		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_distance(pt, ip);
#ifdef GEODEBUG
		printf("dist_ps- distance is %f to intersection point is (%f,%f)\n",
			   *result, ip->x, ip->y);
#endif

		/* otherwise, intersection is not on line segment */
	}
	else
	{
		result = point_distance(pt, &lseg->p[0]);
		tmpdist = point_distance(pt, &lseg->p[1]);
		if (*tmpdist < *result)
			*result = *tmpdist;
		pfree(tmpdist);
	}

	if (ip != NULL)
		pfree(ip);
	pfree(ln);
	return result;
}


/*
 ** Distance from a point to a path
 */
double *
dist_ppath(Point *pt, PATH *path)
{
	double	   *result;
	double	   *tmp;
	int			i;
	LSEG		lseg;

	switch (path->npts)
	{
			/* no points in path? then result is undefined... */
		case 0:
			result = NULL;
			break;
			/* one point in path? then get distance between two points... */
		case 1:
			result = point_distance(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.
			 */
			result = palloc(sizeof(double));
			for (i = 0; i < path->npts - 1; i++)
			{
				statlseg_construct(&lseg, &path->p[i], &path->p[i + 1]);
				tmp = dist_ps(pt, &lseg);
				if (i == 0 || *tmp < *result)
					*result = *tmp;
				pfree(tmp);
			}
			break;
	}
	return result;
}

double *
dist_pb(Point *pt, BOX *box)
{
	Point	   *tmp;
	double	   *result;

	tmp = close_pb(pt, box);
	result = point_distance(tmp, pt);
	pfree(tmp);

	return result;
}


double *
dist_sl(LSEG *lseg, LINE *line)
{
	double	   *result,
			   *d2;

	if (inter_sl(lseg, line))
	{
		result = palloc(sizeof(double));
		*result = 0.0;

	}
	else
	{
		result = dist_pl(&lseg->p[0], line);
		d2 = dist_pl(&lseg->p[1], line);
		if (*d2 > *result)
		{
			pfree(result);
			result = d2;
		}
		else
			pfree(d2);
	}

	return result;
}


double *
dist_sb(LSEG *lseg, BOX *box)
{
	Point	   *tmp;
	double	   *result;

	tmp = close_sb(lseg, box);
	if (tmp == NULL)
	{
		result = palloc(sizeof(double));
		*result = 0.0;
	}
	else
	{
		result = dist_pb(tmp, box);
		pfree(tmp);
	}

	return result;
}


double *
dist_lb(LINE *line, BOX *box)
{
	Point	   *tmp;
	double	   *result;

	tmp = close_lb(line, box);
	if (tmp == NULL)
	{
		result = palloc(sizeof(double));
		*result = 0.0;
	}
	else
	{
		result = dist_pb(tmp, box);
		pfree(tmp);
	}

	return result;
}


double *
dist_cpoly(CIRCLE *circle, POLYGON *poly)
{
	double	   *result;
	int			i;
	double	   *d;
	LSEG		seg;

	if (!PointerIsValid(circle) || !PointerIsValid(poly))
		elog(ERROR, "Invalid (null) input for distance", NULL);

	if (point_inside(&(circle->center), poly->npts, poly->p))
	{
#ifdef GEODEBUG
		printf("dist_cpoly- center inside of polygon\n");
#endif
		result = palloc(sizeof(double));

		*result = 0;
		return result;
	}

	/* 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(&(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(&(circle->center), &seg);
#ifdef GEODEBUG
		printf("dist_cpoly- segment %d distance is %f\n", (i + 1), *d);
#endif
		if (*d < *result)
			*result = *d;
		pfree(d);
	}

	*result -= circle->radius;
	if (*result < 0)
		*result = 0;

	return result;
}	/* dist_cpoly() */


/*---------------------------------------------------------------------
 *		interpt_
 *				Intersection point of objects.
 *				We choose to ignore the "point" of intersection between
 *				  lines and boxes, since there are typically two.
 *-------------------------------------------------------------------*/

static Point *
interpt_sl(LSEG *lseg, LINE *line)
{
	LINE	   *tmp;
	Point	   *p;

	tmp = line_construct_pp(&lseg->p[0], &lseg->p[1]);
	p = line_interpt(tmp, line);
#ifdef GEODEBUG
	printf("interpt_sl- segment is (%.*g %.*g) (%.*g %.*g)\n",
		   digits8, lseg->p[0].x, digits8, lseg->p[0].y, digits8, lseg->p[1].x, digits8, lseg->p[1].y);
	printf("interpt_sl- segment becomes line A=%.*g B=%.*g C=%.*g\n",
		   digits8, tmp->A, digits8, tmp->B, digits8, tmp->C);
#endif
	if (PointerIsValid(p))
	{
#ifdef GEODEBUG
		printf("interpt_sl- intersection point is (%.*g %.*g)\n", digits8, p->x, digits8, p->y);
#endif
		if (on_ps(p, lseg))
		{
#ifdef GEODEBUG
			printf("interpt_sl- intersection point is on segment\n");
#endif

		}
		else
		{
			pfree(p);
			p = NULL;
		}
	}

	pfree(tmp);
	return p;
}


/*---------------------------------------------------------------------
 *		close_
 *				Point of closest proximity between objects.
 *-------------------------------------------------------------------*/

/* close_pl -
 *		The intersection point of a perpendicular of the line
 *		through the point.
 */
Point *
close_pl(Point *pt, LINE *line)
{
	Point	   *result;
	LINE	   *tmp;
	double		invm;

	result = palloc(sizeof(Point));
#ifdef NOT_USED
	if (FPeq(line->A, -1.0) && FPzero(line->B))
	{							/* vertical */
	}
#endif
	if (line_vertical(line))
	{
		result->x = line->C;
		result->y = pt->y;
		return result;

#ifdef NOT_USED
	}
	else if (FPzero(line->m))
	{							/* horizontal */
#endif
	}
	else if (line_horizontal(line))
	{
		result->x = pt->x;
		result->y = line->C;
		return 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(tmp, line);
	return result;
}	/* close_pl() */


/* 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
 */
Point *
close_ps(Point *pt, LSEG *lseg)
{
	Point	   *result;
	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

	result = NULL;
	xh = lseg->p[0].x < lseg->p[1].x;
	yh = lseg->p[0].y < lseg->p[1].y;
	/* !xh (or !yh) is the index of lower x( or y) end point of lseg */

	/* vertical segment? */
	if (lseg_vertical(lseg))
	{
#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)
			return result;

		/* point lines along (to left or right) of the vertical lseg. */

		result = palloc(sizeof(*result));
		result->x = lseg->p[0].x;
		result->y = pt->y;
		return result;
	}
	else if (lseg_horizontal(lseg))
	{
#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)
			return result;

		/* point lines along (at top or below) the horiz. lseg. */
		result = palloc(sizeof(*result));
		result->x = pt->x;
		result->y = lseg->p[0].y;
		return result;
	}

	/*
	 * vert. and horiz. cases are down, now check if the closest point is
	 * one of the end points or someplace on the lseg.
	 */

	/* TODO: Ask if "tmp" should be freed to prevent memory leak */

	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 */
/*	  fprintf(stderr,"below: tmp A %f  B %f   C %f    m %f\n",tmp->A,tmp->B,tmp->C, tmp->m); */
		return 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 */
/*	  fprintf(stderr,"above: tmp A %f  B %f   C %f    m %f\n",tmp->A,tmp->B,tmp->C, tmp->m); */
		return 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);
/*	fprintf(stderr,"tmp A %f  B %f   C %f    m %f\n",tmp->A,tmp->B,tmp->C, tmp->m); */
	result = interpt_sl(lseg, tmp);
/*	fprintf(stderr,"result.x %f  result.y %f\n", result->x, result->y); */
	return result;
}	/* close_ps() */


/* close_lseg()
 * Closest point to l1 on l2.
 */
Point *
close_lseg(LSEG *l1, LSEG *l2)
{
	Point	   *result = NULL;
	Point		point;
	double		dist;
	double	   *d;

	d = dist_ps(&l1->p[0], l2);
	dist = *d;
	memcpy(&point, &l1->p[0], sizeof(point));
	pfree(d);

	if (*(d = dist_ps(&l1->p[1], l2)) < dist)
	{
		dist = *d;
		memcpy(&point, &l1->p[1], sizeof(point));
	}
	pfree(d);

	if (*(d = dist_ps(&l2->p[0], l1)) < dist)
	{
		result = close_ps(&l2->p[0], l1);
		memcpy(&point, result, sizeof(point));
		pfree(result);
		result = close_ps(&point, l2);
	}
	pfree(d);

	if (*(d = dist_ps(&l2->p[1], l1)) < dist)
	{
		if (result != NULL)
			pfree(result);

		result = close_ps(&l2->p[1], l1);
		memcpy(&point, result, sizeof(point));
		pfree(result);
		result = close_ps(&point, l2);
	}
	pfree(d);

	if (result == NULL)
	{
		result = palloc(sizeof(*result));
		memcpy(result, &point, sizeof(*result));
	}

	return result;
}	/* close_lseg() */

/* close_pb()
 * Closest point on or in box to specified point.
 */
Point *
close_pb(Point *pt, BOX *box)
{