<<-- COGHomo.mesa: Basic geometrical operations in homogeneous coordinates>> <<-- last modified by Stolfi - October 15, 1982 4:04 pm>> DIRECTORY COGRandom USING [Toss], COGCart USING [Point, Vector, Box, ScaleFactors]; COGHomo: CEDAR DEFINITIONS IMPORTS COGRandom = BEGIN OPEN Cart: COGCart, Rand: COGRandom; <<-- BASIC DATA TYPES - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - >> Point: TYPE = RECORD [x, y, w: REAL]; Vector: TYPE = Point; << -- The Cartesian coordinates of a point p are [p.x/p.w, p.y/p.w]. The point [0, 0, 0] is indeterminate; if x#0 or y#0, the point [x, y, 0] is at infinity (in the direction of the vector [x, y]). Practically all procedures in this interface return [0, 0, 0] if one or more of the parameters are [0, 0, 0].>> Line: TYPE = Point; << -- The homogeneous line [a, b, c] corresponds to the line a*x+b*y+c = 0 in Cartesian coordinates. An homogeneous point p=[x, y, w] will be on the line if and only if a*x+b*y+c*w=0. A line with a=b=0, c#0 is said to be at infinity (only points at infinity belong to it); if a=b=c=0, the line is indeterminate (any point is in it). A line with c=0 and a#0 or b#0 is a line through the origin. The (homogeneous) vector [a, b, c] is perpendicular to the line [a, b, c].>> << -- The lines r=[a, b, c] and r'=[d*a, d*b, d*c] are essentially equivalent for any d>0. If d<0, the two lines have the same Cartesian coordinates but are considered to have opposite orientations, so that, for example, the distances from a point to those two lines will have opposite signs. The orientation of line r is by definition such that the (Cartesian) vector [a, b] points to the left of r, i.e. the homogeneous point at infinity [a, b, 0] is at the left of r. The line at infinity [0, 0, c] may be thought as a circle of infinite radius centerd at the origin; by definition, its orientation is counterclockwise iff c>0.>> << -- For symmetry reasons, it is convenient to assign an orientation also to homogeneous points. Accordingly, we say that p=[x, y, w] is oriented counterclockwise iff w is nonnegative. It turns out then that the orientation of the line [x, y, w] as seen from the point [x, y, w] agrees with the orientation of the latter.>> << -- TO THINK: maybe this is too meesy; it may be better to leave the orientation of line at infinity undefined, and similarly that of the origin point. It may be better to ignore the orientation altogether, but then we lose the ability to distinguish the two half-planes, and the point at infinity [x, y, 0] becomes the same as [-x, -y, 0], which doesn't suit many applications.>> << -- The orientation of lines allows one to speak unambiguously about their "left" and "right" half-planes. The orientation of points is not that important per se; however, it influences the orientation of points and lines computed from them, so it affects the orientation of derived lines and the direction of points at infinity. For example, the orientation of the line passing through two given points is affected by the orientation of the two points, as well as by their order. As another example, the circumcenter c of three colinear points is at infinity, in a direction perpendicular to the line s supporting the points; the orientation and the order of the points determine on which side of s the center c will lie.>> <<-- PROCEDURES FOR POINTS AND VECTORS - - - - - - - - - - - - - - - - - - - - - - - - - - >> Coords: PROC [a: Point] RETURNS [Cart.Point] = INLINE {RETURN [[a.x/a.w, a.y/a.w]]}; <<-- Converts to Cartesian coordinates. Generates fault if point is infinite or indeterminate.>> InfPt: PROC [u: Point] RETURNS [R: BOOLEAN] = INLINE <<-- Returns TRUE iff point is at infinity (but FALSE for indeterminate point)>> {R _ ((u.w = 0.0) AND (u.x # 0.0 OR u.y # 0.0))}; FinPt: PROC [u: Point] RETURNS [R: BOOLEAN] = INLINE <<--Returns TRUE iff point is finte.>> {R _ u.w # 0.0}; IndPt: PROC [u: Point] RETURNS [R: BOOLEAN] = INLINE <<-- Returns TRUE iff point is indeterminate.>> {R _ (u.w = 0.0 AND u.x = 0.0 AND u.y = 0.0)}; Add: PROC [a, b: Vector] RETURNS [Vector] = INLINE <<-- Adds two homogeneous vectors (or points). The sum of two infinities is indeterminate>> {RETURN [[a.x*b.w+b.x*a.w, a.y*b.w+b.y*a.w, a.w*b.w]]}; Sub: PROC [a, b: Vector] RETURNS [Vector] = INLINE <<-- Subtracts vector b from vector a. The difference of two infinities is indeterminate>> {RETURN [[a.x*b.w-b.x*a.w, a.y*b.w-b.y*a.w, a.w*b.w]]}; Mul: PROC [a: REAL, u: Vector] RETURNS [w: Vector] = INLINE <<-- Multiplies the vector u by a real number a.>> {RETURN [[u.x*a, u.y*a, u.w]]}; Way: PROC [alpha: REAL, p, q: Point] RETURNS [Point]; <<-- Computes the point alpha the way from u to v. If both points are infinite, and alpha#0 and alpha #1, the result is indeterminate.>> Length: PROC [v: Vector] RETURNS [REAL]; <<-- Length of vector v. If v is infinite, the result is LargestNumber; if v is indeterminate, the result is TrappingNaN (don't ask me what that is supposed to mean).>> LengthSq: PROC [v: Vector] RETURNS [REAL]; <<-- Length of vector v, squared (saves a square root extraction). If v is infinite, the result is LargestNumber; if v is indeterminate, the result is TrappingNaN (don't ask me what that is supposed to mean).>> Dist: PROC [p, q: Point] RETURNS [REAL]; <<-- Distance between points p and q. Returns LargestNumber if one point is infinite, and TrappingNaN if both are infinite or one is indeterminate.>> DistSq: PROC [p, q: Point] RETURNS [REAL] = INLINE <<-- Distance between points p and q, squared (saves a square root extraction). Returns LargestNumber if one point is infinite, and TrappingNaN if both are infinite or on is indeterminate.>> {RETURN [LengthSq[Sub[p,q]]]}; <<-- PROCEDURES FOR LINES - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ->> <<-- Predicates for classifying homogeneous lines>> InfLine: PROC [r: Line] RETURNS [R: BOOLEAN] = INLINE <<-- Returns TRUE iff the line is a line at infinity>> {R _ (r.w # 0 AND r.x = 0 AND r.y = 0)}; FinLine: PROC [r: Line] RETURNS [R: BOOLEAN] = INLINE <<-- Returns TRUE iff the line is finite>> {R _ (r.w # 0) AND (r.x # 0 OR r.y # 0)}; IndLine: PROC [r: Line] RETURNS [R: BOOLEAN] = INLINE <<-- Returns TRUE iff the line is indeterminate>> {R _ (r.w = 0.0 AND r.x = 0.0 AND r.y = 0.0)}; Incident: PROC [p: Point, r: Line] RETURNS [R: BOOLEAN]; <<-- Checks if point belongs to line. Returns TRUE if the line or the point are indeterminate (maybe should return FALSE, but TRUE is easier). Note that Incident[p, r] = Incident[r, p]>> ParallelThrough: PROC [r: Line, p: Point] RETURNS [Line] = INLINE <<-- Computes line passing through p and parallel to line r. If the point is at infinity and on the line, or the line is at infinity, the result is indeterminate. If the point is at infinity but off the line, the result is the line at infinity. The orientation of the line is reversed if p.w<0. >> {RETURN [[r.x*p.w, r.y*p.w, -r.x*p.x-r.y*p.y]]}; PerpendicularThrough: PROC [v: Cart.Vector, p: Point] RETURNS [Line] = INLINE <<-- Computes line passing through p and perpendicular to the Cartesian vector v. If the point is at infinity in a direction perpendicular to v, or v is [0, 0], the result is indeterminate. If the point is at infinity but not perp. to v, the result is the line at infinity. The line is oriented 90 degrees clockwise from v (counterclockwise if p.w<0). >> {RETURN [[v.x*p.w, v.y*p.w, -v.x*p.x-v.y*p.y]]}; DistToLine: PROC [p: Point, r: Line] RETURNS [REAL]; <<-- Computes distance from line r to point p. If the line is at infinity or indeterminate, or the point is at infinity in the same direction as the line, returns TrappingNaN If the point is at infinity but in another direction, returns LargestNumber. The sign of the result tells whether the point is to the left of (+), to the right of (-), or on (0) the line r. It is not affected by the orientation of the point.>> SameSide: PROC [p, q: Point, r: Line] RETURNS [R: BOOLEAN]; <<-- Returns TRUE iff the points p and q are on the same side of the line r (or one of them is on the line).>> Intercept: PROC [r, s: Line] RETURNS [Point]; <<-- Computes the intersection of lines r and s. If both lines are coincident (with any orientations) the result is indeterminate. If one line is at infinity, or the lines are distinct but parallel, the result is at infinity (in the direction of the finite lines).The orientation of the point depends on the orientation and position of the lines: it will be counterclockwise iff the counterclockwise angle from r to s, as seen from the intersection point, is in the range [0, \pi].>> Bisector: PROC [p, q: Point] RETURNS [Line]; <<-- Finds the bisector of points p and q. If both points are at infinity or coincident, the line is indeterminate. If only one point is at infinity, the line is at infinity. The bisector is oriented 90 counterclockwise from q-p (reversed if p.w<0 or q.w<0).>> LineThrough: PROC [p, q: Point] RETURNS [Line]; <<-- Finds the line passing through the points p and q. If both points are coincident, or both are at infinity and in opposite directions, the line is indeterminate. If the two points are at infinity but not opposite, the line is at infinity. If one point is at infinity, it will determine the direction of the line. The line is oriented from p to q (reversed if p.w<0 or q.w<0).>> <<-- PROCEDURES FOR DIRECTIONS AND ANGLES - - - - - - - - - - - - - - - - - - - - - - - - >> <<-- An angle between homogeneous vectors, lines or points can be measured and compared by converting it to a Cartesian vector that makes the same angle with the x-axis. This is done by the procedures below. >> VecDir: PROC [v: Vector] RETURNS [Cart.Vector] = INLINE <<-- Gives the direction of the homogeneous vector v as a Cartesian vector Returns [0, 0] if the point is indeterminate or is the origin.>> {RETURN [IF v.w<0 THEN [-v.x, -v.y] ELSE [v.x, v.y]]}; LineNormal: PROC [r: Line] RETURNS [Cart.Vector] = INLINE <<-- Gives the Cartesian vector v perpendicular to the line r and pointing to its left If the line is at infinity or indeterminate, the result is [0, 0].. >> {RETURN [[r.x, r.y]]}; PtPtDir: PROC [p, q: Point] RETURNS [Cart.Vector] = INLINE <<-- Gives the direction of the segment or ray from p to q as a Cartesian vector Returns [0, 0] if the points are coincident or both infinite. Ignores the orientation of the point.>> {RETURN [VecDir[Sub[q, p]]]}; CCW: PROC [p, q, r: Point] RETURNS [BOOLEAN]; <<-- Returns TRUE iff the angle pqr is strictly convex, i.e. iff from some finite point o the points p, q and r are seen counterclockwise in that order.>> CircumCenter: PROC [p, q, r: Point] RETURNS [Point]; <<-- Computes the circumcenter of the points p, q and r. Assumes they are given in counterclockwise order around the triangle pqr.>> <<-- POINT TRANSFORMATIONS FOR PLOTTING - - - - - - - - - - - - - - - - - - - - - - - - - - >> ScalePoint: PROC [p: Point, sf: Cart.ScaleFactors] RETURNS [ps: Cart.Point]; <<-- Transforms point p according to the scale factors sf. If p is infinite, will map it into a point in the same direction (as seen from the scaled origin) and at a very large distance. If p is indeterminate, will map it as if it were [0, 0].>> ScaleVector: PROC [v: Vector, sf: Cart.ScaleFactors] RETURNS [vs: Cart.Point]; <<-- Transforms vector v according to the scale factors sf (will not displace origin). If p is infinite, will map it into a vector in the same direction (as seen from the scaled origin) and at a very large distance. If v is indeterminate, will map it as if it were [0, 0].>> ScaleSeg: PROC [p, q: Point, sf: Cart.ScaleFactors] RETURNS [ps, qs: Cart.Point]; <<-- Transforms the segment pq according to the scale factors sf. If only one of the endpoints is at infinity, will map it into a suitably far away point, so as to preserve the direction of the segment. If only one endpoint is indeterminate, will map it to the other endpoint. If both p and q are infinite or indeterminate, will map them as if they were [0.0, 0.0, 1.0].>> <<-- MISCELLANEOUS PROCEDURES - - - - - - - - - - - - - - - - - - - - - - - - >> Below: PROC [p, q: Point] RETURNS [R: BOOL]; <<-- Returns TRUE if p is strictly below q. Works if one of them is infinite; even if both are infinite, provided p.y and q.y have different signs.>> ToTheLeftOf: PROC [p, q: Point] RETURNS [R: BOOL]; <<-- Returns TRUE if p is strictly to the left of q. Works if one of them is infinite; even if both are infinite, provided p.x and q.x have different signs.>> Precedes: PROC [p, q: Point] RETURNS [R: BOOL]; <<-- Returns TRUE if p precedes q in lexicographical order, i.e. p.x < q.x or p.x = q.x and p.y < q.y (after scaling by p.w and q.w). Works if one of them is infinite; even if both are infinite, provided p.x and q.x have opposite signs. In any case, the answer is transitively consistent.>> InCircle: PROC [p, q, r, s: Point] RETURNS [R: BOOL]; <<-- Checks whether the point s is inside the circle passing through p, q, and r. More specifically, if all of p, q, r, and s are cocircular (or collinear), then it returns FALSE. If p, q, and r are collinear but s is not, it returns TRUE. Else it returns TRUE iff s is strictly inside the circle determined by p, q, and r. This predicate satisfies the identitities:>> <<-- InCircle[p, q, r, s] = InCircle[r, q, p, s] = InCircle[p, s, q, r].>> <<-- If p, q, r, and s are not cocircular, then we also have>> <<-- InCircle[p, q, r, s] = NOT InCircle[q, r, s, p];>> <<-- The computation is done by intersecting the bisectors of pr and qs and checking which pair is closest to the intersection. Assumes p, q, r, s form a simple quadrilateral..>> EvenPerm: PROC [p, q, r: Point] RETURNS [BOOL] = INLINE <<-- Assumes p, q, and r are collinear and distinct. It returns TRUE if pqr in this order form an even permutation of their linear ordering.>> {pq, qr, rp: NAT; pq _ IF Precedes[p, q] THEN 1 ELSE 0; qr _ IF Precedes[q, r] THEN 1 ELSE 0; rp _ IF Precedes[r, p] THEN 1 ELSE 0; RETURN [pq+qr+rp = 2]}; Random: PROC [box: Cart.Box] RETURNS [v: Vector] = INLINE <<-- Returns a random vector in the given box.>> {RETURN [[Rand.Toss [box.min.x, box.max.x], Rand.Toss [box.min.y, box.max.y], 1.0]]}; END...