Define boolean functions for lseg <, <=, <>, >=, >

Define close_ls(), close_lseg(), lseg_length().
Write real code for close_sb(), close_pb(), inter_sb(), inter_lb().
Repair lseg_perp() which determines if two lsegs are perpendicular.
Repair lseg_dt() distance between two lsegs.
Note: close_sl() is clearly broken  but will repair later
 (calculating point on lseg rather than point on line).

src/backend/utils/adt/geo_ops.c

@@ -7,7 +7,7 @@
  *
  *
  * IDENTIFICATION
- *	  $Header: /cvsroot/pgsql/src/backend/utils/adt/geo_ops.c,v 1.29 1998/01/07 18:46:47 momjian Exp $
+ *	  $Header: /cvsroot/pgsql/src/backend/utils/adt/geo_ops.c,v 1.30 1998/02/03 15:55:58 thomas Exp $
  *
  *-------------------------------------------------------------------------
  */
@@ -765,7 +765,6 @@ box_diagonal(BOX *box)
 	p2.x = box->low.x;
 	p2.y = box->low.y;
 	return (lseg_construct(&p1, &p2));
-
 }
 
 /***********************************************************************
@@ -792,7 +791,9 @@ line_construct_pm(Point *pt, double m)
 	result->B = -1.0;
 	result->C = pt->y - m * pt->x;
 
+#if FALSE
 	result->m = m;
+#endif
 
 	return (result);
 } /* line_construct_pm() */
@@ -812,7 +813,9 @@ line_construct_pp(Point *pt1, Point *pt2)
 #ifdef GEODEBUG
 		printf("line_construct_pp- line is vertical\n");
 #endif
+#if FALSE
 		result->m = DBL_MAX;
+#endif
 
 	}
 	else if (FPeq(pt1->y, pt2->y))
@@ -824,7 +827,9 @@ line_construct_pp(Point *pt1, Point *pt2)
 #ifdef GEODEBUG
 		printf("line_construct_pp- line is horizontal\n");
 #endif
+#if FALSE
 		result->m = 0.0;
+#endif
 
 	}
 	else
@@ -840,7 +845,9 @@ line_construct_pp(Point *pt1, Point *pt2)
 		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
+#if FALSE
 		result->m = result->A;
+#endif
 	}
 	return (result);
 } /* line_construct_pp() */
@@ -1249,7 +1256,7 @@ path_inter(PATH *p1, PATH *p2)
 		b2.low.y = Min(p2->p[i].y, b2.low.y);
 	}
 	if (!box_overlap(&b1, &b2))
-		return (0);
+		return (FALSE);
 
 	/* pairwise check lseg intersections */
 	for (i = 0; i < p1->npts - 1; i++)
@@ -1259,16 +1266,18 @@ path_inter(PATH *p1, PATH *p2)
 			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 (1);
+				return (TRUE);
 		}
 	}
 
 	/* if we dropped through, no two segs intersected */
-	return (0);
-}
+	return (FALSE);
+} /* path_inter() */
 
-/* this essentially does a cartesian product of the lsegs in the
-   two paths, and finds the min distance between any two lsegs */
+/* 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)
 {
@@ -1305,7 +1314,7 @@ path_distance(PATH *p1, PATH *p2)
 		}
 
 	return (min);
-}
+} /* path_distance() */
 
 
 /*----------------------------------------------------------
@@ -1500,6 +1509,10 @@ point_distance(Point *pt1, Point *pt2)
 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));
 }
 
@@ -1557,7 +1570,9 @@ lseg_in(char *str)
 		|| (*s != '\0'))
 		elog(ERROR, "Bad lseg external representation '%s'", str);
 
+#if FALSE
 	lseg->m = point_sl(&lseg->p[0], &lseg->p[1]);
+#endif
 
 	return (lseg);
 } /* lseg_in() */
@@ -1586,7 +1601,9 @@ lseg_construct(Point *pt1, Point *pt2)
 	result->p[1].x = pt2->x;
 	result->p[1].y = pt2->y;
 
+#if FALSE
 	result->m = point_sl(pt1, pt2);
+#endif
 
 	return (result);
 }
@@ -1600,9 +1617,24 @@ statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2)
 	lseg->p[1].x = pt2->x;
 	lseg->p[1].y = pt2->y;
 
+#if FALSE
 	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.
  *---------------------------------------------------------*/
@@ -1641,6 +1673,15 @@ lseg_parallel(LSEG *l1, LSEG *l2)
 				 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)
 {
@@ -1650,11 +1691,15 @@ lseg_perp(LSEG *l1, LSEG *l2)
 	m1 = point_sl(&(l1->p[0]), &(l1->p[1]));
 	m2 = point_sl(&(l2->p[0]), &(l2->p[1]));
 
-	if (!FPzero(m1))
-		return (FPeq(m2 / m1, -1.0));
-	else if (!FPzero(m2))
+#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));
-	return (0);					/* both 0.0 */
 } /* lseg_perp() */
 
 bool
@@ -1677,7 +1722,40 @@ lseg_eq(LSEG *l1, LSEG *l2)
 			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() */
 
 
 /*----------------------------------------------------------
@@ -1699,7 +1777,11 @@ lseg_distance(LSEG *l1, LSEG *l2)
 	return (result);
 }
 
-/* distance between l1, l2 */
+/* 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)
 {
@@ -1715,17 +1797,12 @@ lseg_dt(LSEG *l1, LSEG *l2)
 	d = dist_ps(&l1->p[1], l2);
 	result = Min(result, *d);
 	pfree(d);
-#if FALSE
-/* XXX Why are we checking distances from all endpoints to the other segment?
- * One set of endpoints should be sufficient - tgl 97/07/03
- */
 	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);
-#endif
 
 	return (result);
 } /* lseg_dt() */
@@ -2181,10 +2258,14 @@ close_pl(Point *pt, LINE *line)
 } /* close_pl() */
 
 
-/* close_ps -
+/* 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.
  */
 Point *
 close_ps(Point *pt, LSEG *lseg)
@@ -2206,47 +2287,149 @@ close_ps(Point *pt, LSEG *lseg)
 		result = point_copy(&lseg->p[!yh]);
 	else if (pt->y > lseg->p[yh].y)
 		result = point_copy(&lseg->p[yh]);
-	if (result)
+	if (result != NULL)
 		return (result);
-#if FALSE
-	if (FPeq(lseg->p[0].x, lseg->p[1].x))		/* vertical */
-#endif
+
+	/* vertical segment? */
 	if (lseg_vertical(lseg))
 	{
+#ifdef GEODEBUG
+printf("close_ps- segment is vertical\n");
+#endif
+		result = palloc(sizeof(*result));
 		result->x = lseg->p[0].x;
 		result->y = pt->y;
 		return (result);
-#if FALSE
-		}
-		else if (FPzero(lseg->m))
-		{						/* horizontal */
-#endif
 	}
 	else if (lseg_horizontal(lseg))
 	{
+#ifdef GEODEBUG
+printf("close_ps- segment is horizontal\n");
+#endif
+		result = palloc(sizeof(*result));
 		result->x = pt->x;
 		result->y = lseg->p[0].y;
 		return (result);
 	}
 
-#if FALSE
-	invm = -1.0 / lseg->m;
-#endif
 	invm = -1.0 / point_sl(&(lseg->p[0]), &(lseg->p[1]));
 	tmp = line_construct_pm(pt, invm);
 	result = interpt_sl(lseg, tmp);
 	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)
 {
-	/* think about this one for a while */
-	elog(ERROR, "close_pb not implemented", NULL);
+	LSEG		lseg,
+				seg;
+	Point		point;
+	double		dist,
+			   *d;
 
-	return (NULL);
-}
+	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)
 {
@@ -2257,6 +2440,7 @@ close_sl(LSEG *lseg, LINE *line)
 	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)
@@ -2269,14 +2453,88 @@ close_sl(LSEG *lseg, LINE *line)
 	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)
 {
-	/* think about this one for a while */
-	elog(ERROR, "close_sb not implemented", NULL);
+	Point	   *result;
+	Point	   *pt;
+	Point		point;
 
-	return (NULL);
-}
+	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)
@@ -2374,10 +2632,10 @@ on_ppath(Point *pt, PATH *path)
 			b = point_dt(pt, &path->p[i + 1]);
 			if (FPeq(a + b,
 					 point_dt(&path->p[i], &path->p[i + 1])))
-				return (1);
+				return (TRUE);
 			a = b;
 		}
-		return (0);
+		return (FALSE);
 	}
 
 	return (point_inside(pt, path->npts, path->p));
@@ -2404,7 +2662,7 @@ on_ppath(Point *pt, PATH *path)
 		if (FPeq(yh, yl))		/* horizontal seg? */
 			if (FPge(pt->x, xl) && FPle(pt->x, xh) &&
 				FPeq(pt->y, yh))
-				return (1);		/* pt lies on seg */
+				return (TRUE);		/* pt lies on seg */
 			else
 				continue;		/* skip other hz segs */
 		if (FPlt(yh, pt->y) ||	/* pt is strictly below seg */
@@ -2420,7 +2678,7 @@ on_ppath(Point *pt, PATH *path)
 					 &path->p[NEXT(i)]) +
 			path->p[i].x;
 		if (FPeq(x, pt->x))		/* pt lies on this seg */
-			return (1);
+			return (TRUE);
 
 		/* does the seg actually cross the ray? */
 
@@ -2470,25 +2728,108 @@ inter_sl(LSEG *lseg, LINE *line)
 	if (tmp)
 	{
 		pfree(tmp);
-		return (1);
+		return (TRUE);
 	}
-	return (0);
+	return (FALSE);
 }
 
-/* XXX segment and box should be able to intersect; tgl - 97/01/09 */
-
+/* 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)
 {
-	return (0);
-}
+	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);
 
-/* XXX line and box should be able to intersect; tgl - 97/01/09 */
+	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)
 {
-	return (0);
+	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);
 }
 
 /*------------------------------------------------------------------
@@ -3237,7 +3578,7 @@ int4
 poly_npoints(POLYGON *poly)
 {
 	if (!PointerIsValid(poly))
-		return (0);
+		return (FALSE);
 
 	return (poly->npts);
 } /* poly_npoints() */