-- To do: Move some inlines to Impl

-- To do: Think of a better definition for HWay and Clos[p,q].

-- A point of the Two-Sided Plane (2SP) is represented in homogeneous coordinates by a triplet of real coordinates [x,y | w], not all of them zero, with the convention that two triplest [x,y | w] and [x', y' | w'] represent the same point iff we have x=ax', y=ay', and w=aw' for some positive real a.

-- An (oriented) line of the 2SP is also represented by a triplet <a, b | c> of real coordinates, not all three zero. We say that a point [x,y | w] lies on that line, to its left, or to its right, depending on whether ax+by+cw is zero, positive, or negative. We say that a line passes to the left of a point iff th point lies to the left of the line, and vice-versa.

-- Both points and lines are represented in thisinterface as records with three real components, named x, y, z.

-- The natural correspondence between the 2SP and the oriented Cartesian plane is defined in the following way: to the point [x,y | w] there corresponds the point (x/w,y/w); to the line <a,b | c> there corresponds the oriented line with equation ax+by+c=0 and left side defined by the normal direction vector (a,b)/sqrt(a^2+b^2).

-- The natural image L of a line <a,b | c> passes at distance c/sqrt(a^2+b^2) from the Cartesian origin (0,0). The longitudinal orientation of L (that is, the direction in which we must face so that the left side of L lies to our left) is the vector (b, -a)/sqrt(a^2+b^2). It follows that L is oriented counterclockwise as seen from the origin iff c>0. If c=0, L passes through (0,0). The two lines <a, b | c> and <-a,- | -c> map onto the same Cartesian line with opposite orientations; they are said to be the opposite of each other. Two lines <a, b | c> and <a, b | c'> are said to be parallel, while <a,b | c> and <-a,-b | c'> are antiparallel.

-- This correspondence is not defined for the points with w=0 or for the lines with a=b=0. If we consider only the points [x,y | w] with w>0, the natural correspondence is an isomorphism between the 2SP and the oriented Cartesian plane: [x,y | w] is on or to the left of <a,b | c> if and only if the same relation holds for their Cartesian correspondents. The points with w>0 constitute the top side of the 2SP, also called the "homogeneous plane" or "homoplane".

-- The bottom side of the 2SP, also called the "homogeneous antiplane" of just "antiplane", consists of the points with w<0. The natural correspondence also maps the bottom side one-to-one onto the Cartesian plane, but it falls short of being an isomorphism: a point [x,y | w] on the bottom side lies to the left of a line <a,b | c> if and only if the corresponding point lies to the right of the corresponding line.

-- The points [x,y | w] and [-x,-y | -w] of the 2SP that map onto the same Cartesian point (x/w, y/w) are said to be the antipodes of each other. A line that passes through a point also passes through its antipode. The top origin of the 2SP is the point [0,0 | 1] (and its positive multiples); its antipode [0,0 | -1] is called the bottom origin. Under the natural correspondence, both are mapped to the Cartesian origin (0,0). A line <a, b | c> passes through the top (and bottom) origin iff c=0.

-- We can imagine that a point [x,y | 0] of the 2SP corresponds to a point at infinite distance from the origin in the in the direction (x,y). Note that [x,y | 0] is distinct from [-x,-y | 0]. That point lies on all lines of the form <-y, x, c> for all c. The lines <0,0 | 1> and <0,0 | -1> are the two lines at infinity of the 2SP, and consist of exactly all the points at infinity. The left side of <0,0 | 1> is the entire top side of the 2SP, whereas that of <0,0 | -1> is the bottom side.

-- The triplets [0, 0 | 0] and <0,0 | 0> are called the indeterminate point and indeterminate line, respectively. They are not considered to be true elements of the 2SP, but they are returned by many operations in this interface when the latter are applied to invalid combinations of arguments. For example, the procedure that computes the intersection of two lines returns [0,0 | 0] when the lines are coincident. Many operations will return [0,0 | 0] or <0, 0 | 0> if one or more of the arguments is [0,0 | 0] of <0,0 | 0>. However, in some cases an arithmetic trap (divide by zero, or such) will result instead, so be careful.

-- An oriented line defines a cyclical ordering of the points lying on it: given any three points p1,q,p2, on L, we can tell whether q lies on the way from p1 to p2 along L (HOW?). This establishes an homeomorphism between any oriented line of the 2SP and an oriented circle. The set of all q's on L on the way from p1 to p2 constitutes by definition the arc on L from p1 to p2. This arc is proper if no two of its points (including the extremeities) are antipodes of each other. Let p1, p2 be two distinct, non-antipodal points of th 2SP; there are exactly two lines L1, L2 passing through those two points, and they are opposite to each other. We define the oriented line p1p2 as the one of the two such that the arc from p1 to p2 of is proper, i.e. the one that goes from p1 to p2 by the shortest way. The oriented line p2p1 is the opposite of the line p1p2 (or both are undefined).

-- Two lines L1, L2 of the 2SP that are not coincident or opposite of each other meet at exactly two antipodal points. Those points will be both at infinity if the lines are parallel; if not, they are both finite.

-- The point [x,y | w] and the line <x,y | w> are dual of each other. If the point p lies on (resp. to the left of) the line L, then the point L* that is dual of L lies on (resp. to the left of) the line p* that is dual of p. Two points are antipodal of each other iff theis duals are opposite lines.

-- Therefore, If L1, L2 are not coincident or opposite of each other, their dual points L1* and L2* define a unique oriented line L1*L2* passing through both points. The dual (L1*L2*)* of this line is a point that lies on both L1 and L2, and is by definition the intersection point L1L2. Note that the point L2L1 is the antipode of the point L1L2

-- Two lines L1, L2 of the 2SP that are not coincident or opposite of each other meet at exactly two antipodal points. Those points will be both at infinity if the lines are parallel; if not, they are both finite. We define the intersection p of L1 and L2 (in that order) as the one of those two points such that we can get L2 by rotating L1 counterclockwise around p by less than 180 degrees.

-- These operations and predicates are "homogeneous", in the sense that the origin and the line at infinity play no special role in them.

-- Returns TRUE iff point is indeterminate.

-- Returns TRUE iff line is indeterminate.

-- Returns the bottom equivalent tof a top point, or the top equivalent of a bottom one.

-- Returns the line r oriented the opposite way.

-- Checks if point belongs to line. Returns TRUE if the line or the point are indeterminate. Note that Incident[p, r] = Incident[r*, p*].

-- The optional parameter tol is the rounding error tolerance; TRUE is returned if the unnormalized signed distance is less than 2^tol.

-- Returns TRUE if p is to the left of the line r. Takes into accout the topness or botomness of p. Returns FALSE if the point and/or line are indeterminate. Note that LeftOf[p, r] = LeftOf[r*, p*]

-- Returns TRUE if p is to the right of the line r. Takes into accout the topness or botomness of p. Returns FALSE if the point and/or line are indeterminate. Note that RightOf[p, r] = RightOf[r*, p*]

-- Returns TRUE iff the points p and q are on the same side of the line r, or one of them is on r, or both of them are on r.

-- Returns TRUE iff the points p and q are on the same side of the line r, and neither of them is on r.

-- Finds the line passing through the points p and q. If the two points are coincident or one is indeterminate, the line is indeterminate. Otherwise, the line is oriented in the direction that gives the ``shortest'' path from p to q.

-- Note that PPLine[p,q] = Opp[PPLine[q,p]] = PPLine[Anti[p], Anti[q]] =Opp[PPLine[p, Anti[q]] = Opp[PPLine[Anti[p], q].

-- Computes the intersection of lines r and s. If both lines are coincident or opposite the result is indeterminate. If one line is at infinity, or the lines are distinct but parallel or antiparallel, the result is at infinity.

-- The intersection of two lines, if defined, consists of two antipodal points. The returned point will be the one such that s can be obtained by turning r counterclockwise around it by less than 180 degrees.

-- Note that Intercept[r,s] = PPLine[r*, s*]. Note also that Intercept[r,s] = Anti[Intercept[s,r]] = Intercept[Opp[r], Opp[s]] =Anti[Intercept[r, Opp[s]] = Anti[Intercept[Opp[r], s].

-- Returns TRUE iff the triangle pqr (where each side is a proper line segment) is striclty counterclockwise oriented. We have PosOr[p, q, r] = LeftOf [p, LineThru[q,r]] = LeftOf [q, LineThru[r,p]] = LeftOf [r, LineThru[p,q]].

-- Returns TRUE iff the triangle pqr (where each side is a proper line segment) is striclty clockwise oriented. We have NegOr[p, q, r] = RightOf [p, PPLine[q,r]] = RightOf [q, PPLine[r,p]] = RightOf [r, PPLine[p,q]].

-- Computes the point alpha the way from u to v, in the homogeneous way. The result lies on the line uv, but the Cartesian diatances from the two points are not proportional to alpha and 1-alpha. In particular, if u and v are on opposite sides of the plane, as alpha goes from 0 to 1 the interpolated point will start at u, move faster and faster away from v and off to infinity, and then it will reappear on the other side of the plane, moving still on the line uv from minus infinity towards v.

-- In particular, if u and v are antipodes of each other, will return u if alpha<1/2, the indeterminate point if alpha=1/2, and v if alpha >1/2. If u is on top and v = -u, will give the origin for alpha=1/2. If u is finite and v = -Anti[u], will give a point at infinity for alpha=1/2.

-- Computes the ``homogeneous'' closeness between p and q. HClos[p,q] is an alteernate non-euclidean measure of the distance between p and q, that is better behaved than Dist and Clos.

-- If one of the points is indeterminate, an hardware error (divide by zero) occurs. HClos[p,q] is 1 when q=p, is positive and decreasing as q moves away from p towards the line p*, is zero when q is on p*, is negative and decreasing as q moves on the other side of p* towards Anti[p], and is -1 when q=Anti[p]. Similarly, if q is finite the same statements hold with p and q interchanged.

-- The set S(p,c) of all points q satisfying HClos[p,q]=c is a cubic (ellipse, parabola or hiperbole), not a circle. Therefore, HClos must not be used to compare the distances of two points q1 and q2 from p, unless q1, q2, and p lie on the same line. However, each sets S(p,c) is topologically equivalent to a circle separating $p$ from Anti[p], and as c goes from 1 to -1 the set S(p,c) expands away from p and towards its antipode. Therefore HClos[p,q] can be used as a metric in the 2SP. Note that it is more defined than Clos.

-- The set S(p,c) of all points q satisfying HClos[p,q]>c is a strictly convex subset of the 2SP iff c>0, and strictly concave otherwise. Indeed, S(p,c) is the antipode of the complement of S(p, -c).

-- Homogeneous closeness of p and q, squared (saves a square root extraction). If one of the points is indeterminate, an hardware error (divide by zero) occurs. The sign is always positive or zero. HClosSq decreases from 1 to 0 as q moves from p to the line p*, and then back to 1 as correctly orders the points q on the same side as p according to increasing distance from p, including the points at infinity. Points q on the other side have the same ClosSq as their antipodes. If q is finite, the same statements hold with p and q interchanged.

-- Computes the ``homogeneous'' distance from line r to point p, which is the same as the homogeneous closeness between r* and p. If one of the arguments is indeterminate, an hardware error (divide by zero) occurs. Otherwise HLPDist[r,p] decreases from +1 when p=r* to 0 when p is on r to -1 when p goes to Anti[r*]. The sign of the result is positive, negative, or zero, depending on whether the point is to the left of r, to the right of r, or on r, respectively. In particular, HLPDist[r, p]= - HLPDist[r, Anti[p]].

-- Computes the ``homogeneous'' distance from line r to point p, squared (saves a square root extraction). If one of the arguments is indeterminate, an hardware error (divide by zero) occurs. Otherwise HLPDistSq[r,p] decreases from +1 when p=r* to 0 when p is on r, and increases again to +1 when p goes to Anti[r*]. The result is independent of the orientation of r and p.

-- The operations below treat the origin and the lines at infinity in a special way.

-- POINTS

-- Returns TRUE iff point is finte.

{RETURN[u.w # 0.0]};

-- Returns TRUE iff point is at infinity (but FALSE for indeterminate point)

{RETURN[(u.w = 0.0) AND (u.x # 0.0 OR u.y # 0.0)]};

-- Returns TRUE if p is strictly on the top side, i.e. iff p.w>0.

{RETURN [p.w > 0.0]};

-- Returns TRUE if p is strictly on the bottom side, i.e. iff p.w<0.

{RETURN [p.w < 0.0]};

-- Converts homogeneous to Cartesian coordinates. A point p and its antipode are mapped to the same cartesian point. If the point is at infinity or indeterminate, an hardware error (divide by zero) will result.

-- Converts Cartesian to homogeneous coordinates. The result is always on the top plane (w = 1).

{RETURN [[p.x, p.y, 1.0]]};

-- Computes the point alpha the way from u to v, that is, a point on the line pq whose Cartesian distances from p and q are in the ratio alpha:(1-alpha). This actually specifies two points, one on top and one on the bottom side of the 2SP. The result will be on the top side is p and q are on both on the top side or both on the bottom side. The result will be on the bottom side if one of the points is on the top and the other is on the bottom.

-- In particular, for points u,v on the top side, if alpha=0 the result is u, if alpha=1 the result is v, if alpha = 1/2 the result is the midpoint of the segment uv. As alpha goes to infinity, the result will also go to infinity. If the two points are finite but coincident, the result is independent of alpha. If q is infinite but p is finite, Way[p, q, alpha] will be p if alpha=0, q if alpha>0, and -q if alpha<0. If one point If both points are infinite, or one is indeterminate, the result is indeterminate. Note that Way[p, q, alpha] = Way [q, p, 1-alpha].

-- Cartesian distance between points p and q. If one of the points is infinite, or indeterminate, an hardware error (divide by zero) occurs. Note: Dist[p, q] is zero if p=q or p=Anti[q]; otherwise it is positive or negative depending on whether p and q are on the same or opposite sides of the 2SP.

-- Cartesian distance between points p and q, squared (saves a square root extraction). If one of the points is infinite or indeterminate, an hardware error (divide by zero) occurs. The sign is always positive or zero.

-- Closeness of p and q. Roughly speaking, Clos[p,q] is to Dist[p,q] what Shrt[v] is to Length[v]. It gives a non-Euclidean measure of how close the two points are to each other, and is better behaved than Dist[p,q] outside the top plane 2SP.

-- If one of the points is indeterminate, or both are infinite, an hardware error (divide by zero) occurs. If p is finite, Clos[p,q] is 1 when q=p, is positive and decreasing as q moves away from p on the same side of the 2SP, is zero when q is at infinity, is negative and decreasing as q moves on the other side towards Anti[p], and is -1 when q=Anti[p]. The points q satisfying Clos[p,q]=c constitute a circle of radius r=SqRt[1-c^2]/c centered at p (a negative radius means the circle is on the side opposite to p). Similarly, if q is finite the same statements hold with p and q interchanged.

-- Therefore, Clos[p,q] correctly sorts the points q of the 2SP according to their ``topological'' distance from p, and vice-versa.

-- Closeness of p and q, squared (saves a square root extraction). If one of the points is indeterminate, or both are infinite, an hardware error (divide by zero) occurs. The sign is always positive or zero. If p is finite, ClosSq correctly orders the points q on the same side as p according to increasing distance from p, including the points at infinity. Points q on the other side have the same ClosSq as their antipodes. If q is finite, the same statements hold with p and q interchanged.

-- VECTORS

-- A vector on the 2SP is formally the same as a point, except that (for documentation purposes) addition and subtraction are defined only for vector-vector, point-vector, and vector-point pairs. Also, the length function and the product by a scalar are defnedfor vectors, not points. We say that the point [x,y | w] is the tip of the vector [x,y | w].

-- Two vectors [x,y | w] and [x', y' | w'] are said to be parallel if x'=ax and y'=ay for some positive constant a, and antiparallel if that is true for some negative constant a. The opposite of a vector v = [x,y | w] is by definition -v = [-x,-y | w] Note that v and Anti[v] = [-x, -y | -w] are antiparallel (but not oposite) even though their tips map to the same Cartesian point.

-- Adds two homogeneous vectors (or points). The sum of two infinities is indeterminate. Topness and bottomness act as + and - under multiplication: the sum of two top or two bottom vectors is a top vector, the sum of a top and a bottom is a bottom.

-- Subtracts vector b from vector a. The difference of two infinities is indeterminate. Top and bottom as in Add.

-- Multiplies the vector u by a real number a. Topness not affected.

-- Cartesian length of vector v. If v is infinite, or indeterminate, an hardware error (divide by zero) occurs. Note: the sign of Length[v] is negative if v is nonzero and on the bottom side; it is zero if v is the origin or the antiorigin.

-- Cartesian length of vector v, squared (saves a square root extraction). If v is infinite or indeterminate, an hardware error (divide by zero) occurs. The sign is always positive or zero.

-- If v is the origin, the anti-origin, or indeterminate, returns the indeterminate vector. Otherwise returns a unit vector on the top side, pointing from the origin towards v by the shortest route. That is, if v is on the top side or at infinity, returns a unit vector parallel to v; if v is on the bottom side, returns a unit vector antiparallel to Anti[v].

-- Shortness of v, a measure of how close v is to the origin. Basically SqRt[1/(1+Length[v]^2)], except fo the sign and range of definition. If v is indeterminate, an hardware error (divide by zero) occurs. Otherwise the shortness of v decreases monotonically from 1 to 0 as the length of v goes from zero to infinity on the top plane, and continues decreasing from 0 to -1 as v goes from infinity towards the antiorigin on the bottom side. Therefore, Shrt correctly orders the points of the whole 2SP by "topological" distance from the origin. The equation Shrt[v]=c, for 1>c>-1, defines a circle centered at the origin with radius r=SqRt[1-c^2]/c (a negative radius means the circle is on the antiplane).

-- Shortness of vector v, squared (saves a square root extraction). If v is indeterminate, an hardware error (divide by zero) occurs. The sign is always positive or zero. ShrtSq correctly orders the points of the top side of the 2SP according to increasing distance from the origin, including the points at infinity. Points on the bottom side have the same ShrtSq as their antipodes.

-- LINES

-- Returns TRUE iff line is finte.

{RETURN[r.x # 0.0 OR r.y # 0.0]};

-- Returns TRUE iff line is at infinity (but FALSE for indeterminate line)

{RETURN[r.w # 0.0 AND r.x = 0.0 AND r.y = 0.0]};

-- Computes the signed distance from line r to point p. If the line is at infinity or indeterminate, or the point is at infinity, an hardware error (divide by zero) occurs. The sign of the result tells whether the point is to the left of (+), to the right of (-), or on (0) the line r. Takes into account the side in which the point lies, so LPDist[r, p]= - LPDist[r, Anti[p]].

-- Finds the bisector of points p and q. If either point is indeterminate, or both points are at infinity, or the points are coincident, the line is indeterminate. If only one point is at infinity, the line is at infinity. The bisector is oriented 90 counterclockwise from the line pq as seen from the midpoint of the proper segment pq.

-- LINES AND VECTORS

-- The normal vector of a line <a, b | c> is by definition [a, b | SqRt[a^2+b^2]]. This vector is indeterminate if a=b=0, that is, the line <a,b | c> is at infinity or indeterminate. The notmal vector is perpendicular to the line. If we add the normal vector to a point on the top half of the line, we will obtain a point still on the top plane and to the left of the line.

-- If the line is at infinity or indeterminate, returns the zero (top) vector. Otherwise returns an unnormalized vector on the top side, directed 90 counterclockwise from the line r (as seen from the top plane).

-- If the line is at infinity or indeterminate, returns the zero (top) vector. Otherwise returns an unnormalized vector on the top side parallel to r (as seen from the top plane).

-- Computes the line with normal vector parallel to v and passing through p. If v is zero, or p is indeterminate, or p is at infinity in a direction orthogonal to v, the result is indeterminate. If p is at infinity but in some other direction, the result is one of the lines at infinity. The vector v does not have to be normalized.

-- 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.

-- The procedures below compare the coordinates of the cartesian images of two points. They work properly only if both points are finite; if not, an hardware error (divide by zero) will result.

-- Returns TRUE if the y coordinate of the Cartesian image of p is greater than that of q.

-- 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.

-- 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.

-- 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]]};

-- 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]]};

-- 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]]]};

-- Computes the circumcenter of the points p, q and r. Assumes they are given in counterclockwise order around the triangle pqr.

-- 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].

-- 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].

-- 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].

-- 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..

-- 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]};

-- 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]]};