-- COGJunk.mesa: Old junk from Voronoi/Delaunay Diagrams
-- last modified by Stolfi September 10, 1983 10:04 pm
DIRECTORY
Rope USING [ROPE];
COGJunk: CEDAR DEFINITIONS =
BEGIN
-- PRIMITIVES FOR DIAGRAM CONSTRUCTION - - - - - - - - - - - - - - - - - - - - - - - - - - -
VertexCount, FaceCount: NAT ← 0;
MakeVertex: ENTRY PROC [anEdge: Edge ← NIL] RETURNS [v: Vertex] =
BEGIN
v ← NEW[VertexRec ←
[anEdge: anEdge,
pos: [0.0, 0.0, 0.0],
no: (VertexCount ← VertexCount+1),
data: NIL]]
END;
MakeFace: ENTRY PROC [anEdge: Edge ← NIL] RETURNS [f: Face] =
BEGIN
f ← NEW[FaceRec ←
[anEdge: anEdge,
pos: [0.0, 0.0, 0.0],
no: (FaceCount ← FaceCount+1),
data: NIL]]
END;
ShowDualOfMe: Menus.MenuProc =
BEGIN
-- parameters: [viewer: Viewer, clientData: REF ANY, redButton: BOOL]
 vData: ViewerData ← NARROW [viewer.data];
out.PutF["Called ShowDualOfMe...\n"];
vData.showDelaunay ← NOT vData.showDelaunay;
vData.action ← $Repaint;
TellMother[]
END;
ComputeVoronoiVertex: INTERNAL PROC [p, q, r: Homo.Point] RETURNS [v: Homo.Point] = TRUSTED
BEGIN
-- Computes the Voronoi vertex that corresponds to the Delaunay face passing through p, q and r (i.e., their circumcenter). If p, q and r do not form a convex angle, they are assumed to be on the convex hull; in that case, the infinity point perpendicular to pq is returned.
 IF Homo.Convex [p, q, r]
THEN {v ← Homo.CircumCenter [p, q, r]}
ELSE {dir: Homo.Vector ← Homo.Sub[q, p];
v ← IF dir.w < 0.0 THEN [dir.y, -dir.x, 0.0] ELSE [-dir.y, dir.x, 0.0]}
END;
ConnectPoints: INTERNAL PROC [p, q: REF DVertex, ep, eq: REF DEdge, lv, rv: REF VVertex]
RETURNS [e: REF DEdge] = TRUSTED
BEGIN
-- inserts an edge e (and its symmetrical) between the Delaunay vertices p and q.
 -- ep must be the edge out of p preceding immediately e c.c.w.,
 -- and eq must be the edge out of q immediately preceding e.sym c.c.w.
 out.PutF["Adding edge (%g, %g)", int[p.no], int[q.no]];
IF ep # NIL THEN out.PutF[" after (%g, %g) at p ", int[Org[ep].no], int[Dest[ep].no]];
IF eq # NIL THEN out.PutF[" after (%g, %g) at q", int[Org[eq].no], int[Dest[eq].no]];
out.PutF["...\n"];
e ← NEW [DEdge ← [destV: q, fNext: NIL, vNext: NIL, sym: NIL, leftVV: lv]];
e.sym ← NEW [DEdge ← [destV: p, fNext: NIL, vNext: NIL, sym: e, leftVV: rv]];
lv.anEdge ← e;
rv.anEdge ← e.sym;
IF ep # NIL THEN
{e.vNext ← ep.vNext; e.sym.fNext ← ep;
ep.vNext.sym.fNext ← e; ep.vNext ← e}
ELSE
{e.vNext ← e; e.sym.fNext ← e; p.anEdge ← e};
IF eq # NIL THEN
{e.sym.vNext ← eq.vNext; e.fNext ← eq;
eq.vNext.sym.fNext ← e.sym; eq.vNext ← e.sym}
ELSE
{e.sym.vNext ← e.sym; e.fNext ← e.sym; q.anEdge ← e.sym}
END;
DeleteEdge: INTERNAL PROC [e: REF DEdge] = TRUSTED
BEGIN
-- Deletes the edge e (and its symmetrical).
 p: REF DVertex = Org[e];
q: REF DVertex = Dest[e];
ep: REF DEdge = VPrev[e];
eq: REF DEdge = VPrev[e.sym];
out.PutF["Deleting edge (%g, %g)...\n", int[p.no], int[q.no]];
IF ep # e THEN
{ep.vNext ← e.vNext; e.vNext.sym.fNext ← ep;
IF p.anEdge = e THEN p.anEdge ← ep}
ELSE
p.anEdge ← NIL;
IF eq # e THEN
{eq.vNext ← e.sym.vNext; e.sym.vNext.sym.fNext ← eq;
IF q.anEdge = e.sym THEN q.anEdge ← eq}
ELSE
q.anEdge ← NIL
END;
DumbLocateInVoronoi: INTERNAL PROC [p: Point, Diag: REF DelaunayDiagram]
RETURNS [np: REF DVertex] = TRUSTED
BEGIN
-- Returns the face of the Voronoi containing point p
 -- Temporary version - uses exhaustive enumeration
 pts: LIST OF REF DVertex ← Diag.points;
pa: REF DVertex;
hp: HPoint ← [p.x, p.y, 1.0];
minD: REAL;
paD: REAL;
np ← pts.first; pts ← pts.rest;
minD ← Homo.DistSq[np.dp, hp];
WHILE pts # NIL DO
pa ← pts.first;
IF Homo.FinPt[pa.dp] AND minD > (paD ← Homo.DistSq[pa.dp, hp]) THEN
{minD ← paD; np ← pa};
pts ← pts.rest
ENDLOOP
END;
ShowEdge: PROC
[Diag: REF DelaunayDiagram, e: REF DEdge, style: GraphicStyle] = TRUSTED
BEGIN
vData: ViewerData = NARROW [Diag.viewer.data];
delau:BOOL = vData.showDelaunay;
p: REF DVertex = Org[e];
q: REF DVertex = Dest[e];
thing: REF GraphicThing;
IF (IF delau THEN Homo.FinPt[p.dp] OR Homo.FinPt [q.dp]
ELSE Homo.FinPt[p.dp] AND Homo.FinPt [q.dp]) THEN
{thing ← NEW [GraphicThing ←
[style: style,
obj: NEW [LineSegment ←
IF delau THEN [org: Org[e].dp, destV: Dest[e].dp]
ELSE [org: LeftVV [e].vp, destV: RightVV [e].vp]]
]];
vData.extraThings ← CONS [thing, vData.extraThings];
ViewerOps.PaintViewer [Diag.viewer, client, FALSE, thing];
Rest[]}
END;
ShowRegion: INTERNAL PROC
[Diag: REF DelaunayDiagram, p: REF DVertex, style: GraphicStyle] = TRUSTED
BEGIN
vData: ViewerData = NARROW [Diag.viewer.data];
delau:BOOL = vData.showDelaunay;
e0, ea: REF DEdge;
prevInf: BOOL ← TRUE;
vs: Polygon ← NIL;
thing: REF GraphicThing;
IF NOT Homo.FinPt[p.dp] THEN RETURN;
thing ← NEW [GraphicThing ← [style: style, obj: NIL]];
e0 ← ea ← p.anEdge;
DO
IF Homo.FinPt [Dest[ea].dp] THEN
{IF prevInf THEN vs ← CONS [RightVV[ea].vp, vs];
vs ← CONS [LeftVV[ea].vp, vs];
prevInf ← FALSE}
ELSE
{prevInf ← TRUE};
ea ← VNext[ea]; IF ea = e0 THEN EXIT
ENDLOOP;
thing.obj ← vs;
vData.extraThings ← CONS [thing, vData.extraThings];
ViewerOps.PaintViewer [Diag.viewer, client, FALSE, thing];
Rest[]
END;
ShowPoint: INTERNAL PROC
[Diag: REF DelaunayDiagram, pt: Homo.Point, style: GraphicStyle] = TRUSTED
BEGIN
vData: ViewerData = NARROW [Diag.viewer.data];
thing: REF GraphicThing =
NEW [GraphicThing ← [style: style, obj: NEW [Homo.Point ← pt]]];
vData.extraThings ← CONS [thing, vData.extraThings];
ViewerOps.PaintViewer [Diag.viewer, client, FALSE, thing];
Rest[]
END;
MakeDelaunay: ENTRY PROC RETURNS [Diag: REF DelaunayDiagram] = TRUSTED
BEGIN
?? rewrite!
-- returns a new Delaunay diagram with only three dummy vData points
 -- at the Mercedes infinities.
 sqrt3: REAL = SqRt[3.0];
pi1: REF DVertex = NEW [DVertex ← [dp: [0.0, sqrt3/3.0, 0.0], no: -1, anEdge: NIL]];
pi2: REF DVertex = NEW [DVertex ← [dp: [-0.5, -sqrt3/6.0, 0.0], no: -2, anEdge: NIL]];
pi3: REF DVertex = NEW [DVertex ← [dp: [0.5, -sqrt3/6.0, 0.0], no: -3, anEdge: NIL]];
v123: REF VVertex = NEW [VVertex ← [vp: [0.0, 0.0, 1.0], anEdge: NIL]];
v321: REF VVertex = NEW [VVertex ← [vp: [0.0, 0.0, 0.0], anEdge: NIL]];
e12: REF DEdge = ConnectPoints [pi1, pi2, NIL, NIL, v123, v321];
e23: REF DEdge = ConnectPoints [pi2, pi3, e12.sym, NIL, v123, v321];
e31: REF DEdge = ConnectPoints [pi3, pi1, e23.sym, e12, v123, v321];
nDV ← 0;
RETURN [NEW [DelaunayDiagram ←
[points: LIST [pi1, pi2, pi3],
viewer: NIL]]]
END;
Mother: PROC [viewer: ViewerClasses.Viewer] = TRUSTED
BEGIN
ENABLE ABORTED => GO TO bye;
vData: ViewerData ← NARROW [viewer.data];
WHILE TRUE DO
WaitStatusChanged;
IF NOT vData.live THEN GO TO bye;
vData.extraThings ← NIL;
ViewerOps.PaintViewer [viewer, client, FALSE, NIL];
IF vData.action = $LeftDown THEN
{LockAndInsertInVoronoi[vData.pt, vData.Diag]}
ELSE IF vData.action = $RightDown THEN
{out.PutF["Sorry, can't delete yet.\n"]}
ELSE IF vData.action = $Repaint THEN
{}
ELSE IF vData.action = $Sort THEN
{LockAndSortRegions [vData.Diag]}
ELSE IF vData.action = $Throw THEN
{n: NAT;
out.PutF ["How many points? "];
n ← in.GetInt[]; out.PutF["\n"];
THROUGH [1..n] DO
LockAndInsertInVoronoi [ThrowAPoint [], vData.Diag]
ENDLOOP;
};
vData.action ← NIL
ENDLOOP;

EXITS bye =>
{Close[in]; Close[out]};
END;
Rest: PROC = TRUSTED {Process.Pause[Process.MsecToTicks[300]]};

[in,out] ← CreateTTYStreams["Voronoi.TTY"];
[] ← Random.Init[,-1];

out.PutF["Hello!\n"];
Process.Detach[FORK Mother[MakeViewer[]]]
END.

-- The vertex records contain the (homogeneous) coordinates of a point that is used when drawing the diagram on the screen. These coordinates are used when displaying the diagram on a Drawing viewer; please note that the topology of the diagram, and in particular the ccw orderings of edges, are completely independent of them (even though they will agree in many applications).
-- The face records too contain the coordinates of a point, that is used to represent the face when drawing the dual of the diagram. As in the case of the vertices, there is no restriction on the coordinates of these points; they need not belong to the corresponding faces, and their positions need not be consistent with the ccw orderings and left/right relationships.
-- The vertices and faces are numbered sequentially upon creation, for debugging and to help in the drawing (the drawing routine uses them to avoid painting an edge twice: it paints only those edges with Dest[e].no >= Org[e].no.)
DelaunayPainter: Draw.PainterProc;
???MakeDelaunayObject: PROC [dl: Delaunay, color: Draw.Color, Voronoi: BOOL ← FALSE]
RETURNS [obj: Draw.Object] =
{RETURN [Draw.Make
[DelaunayPainter,
[data: diag, color: color, style: IF Voronoi THEN 1 ELSE 0]]]};
-- Makes a Delaunay diagram (pointed to by the objData) into a Drawing object that will paint itself either in dual or in primal form, depending on the Voronoi flag.

-- inserts an edge e (and its symmetrical) between the Delaunay vertices p and q.
 -- ep must be the edge out of p preceding immediately e c.c.w.,
 -- and eq must be the edge out of q immediately preceding e.sym c.c.w.
 -- out.PutF["Adding edge (%g, %g)", int[p.no], int[q.no]];
 IF ep # NIL THEN out.PutF[" after (%g, %g) at p ", int[Org[ep].no], int[Dest[ep].no]];
IF eq # NIL THEN out.PutF[" after (%g, %g) at q", int[Org[eq].no], int[Dest[eq].no]];
out.PutF["...\n"];
e ← NEW [DEdge ← [destV: q, fNext: NIL, vNext: NIL, sym: NIL, leftVV: lv]];
e.sym ← NEW [DEdge ← [destV: p, fNext: NIL, vNext: NIL, sym: e, leftVV: rv]];
lv.anEdge ← e;
rv.anEdge ← e.sym;
IF ep # NIL THEN
{e.vNext ← ep.vNext; e.sym.fNext ← ep;
ep.vNext.sym.fNext ← e; ep.vNext ← e}
ELSE
{e.vNext ← e; e.sym.fNext ← e; p.anEdge ← e};
IF eq # NIL THEN
{e.sym.vNext ← eq.vNext; e.fNext ← eq;
eq.vNext.sym.fNext ← e.sym; eq.vNext ← e.sym}
ELSE
{e.sym.vNext ← e.sym; e.fNext ← e.sym; q.anEdge ← e.sym}
AddEdge: PROC [Dp, Dq: Delaunay, ep, eq: REF DEdge, curr: Access] RETURNS [e: REF DEdge];
-- Inserts an edge e (and its symmetrical) between the Delaunay vertices p = Dest[ep] and q = Dest[eq]. Will not compute the left and/or right Voronoi vertirces of e!
-- If Dp and/or Dq are currently part of some Drawing, it is the client's responsibility to Acquire write rights to the latter before calling this procedure.
-- The .parms.data field of both Dp and Dq will be set to e (this is necessary in the case when Dp and/or Dq had no edges before). If Dp and Dq are distinct but are in the object list of the same Drawing, the client should remove one of them before relinquishing the write rights; otherwise, the combined diagram may be painted twice on the screen.
DeleteEdge: PROC [De: Delaunay, e: REF DEdge];
-- Deletes the edge e (and its symmetrical) from the given Delaunay. Will not compute the left and/or right Voronoi vertirces of e!
-- The client should have write access rights to the Drawing in which De is being displayed (if it is the case). If the deletion of e causes

I hope it is true that any undirected graph for which such orderings can be defined can be embedded in some two-dimensional manifold in such a way that the orderings and incidences are preserved, no edges intersect, every vertex is contractible to a point, and every face is isomorphic to a disk. And vice-versa.

I also hope that the necessary and sufficient condition for that manifold to be isomorphic to a sphere (or to the plane plus the infinity point) is that the diagram be conneted and that the Euler relation for the sphere be satisfied, namely v-e+f=2. In general, a diagram with c components and v-e+f=2c-2h can be embedded in some 2-D manifold consisting of c disjoint spheres with a total of h handles.
-- The dual of a diagram is obtained by exchanging vertices with faces, and reversing the ccw orderings. This is the same to say that the FNext[e] in the dual diagram will be the same as VPrev[e] in the original, the origin of e in the dual will be its left face, and so forth.
-- A diagram (or rather, a connected component of it) is referred to by giving a pointer to one of its edges. The trivial diagram (one face, one vertex, zero edges) is represneted by a special "dummy" edge, with the sym, fNext and vNext fields pointing to itself, and whose function is only to hold the pointers to the vertex and face records.
-- An edge record also contains a pointer org to its origin vertex, and a pointer left to the face at its left. The representation of vertices and faces is client-defined; the routines in this module treat those pointers as REF ANY.
EdgeRec: TYPE = RECORD -- One half of an undirected edge
 [fNext: Edge,   -- directed edge following e in ccw order as seen from Left[e]
 vNext: Edge,   -- directed edge that follows e in ccw order as seen from Org[e]
sym: Edge,    -- record for the antiparallel edge e'
 org: Vertex ← NIL,  -- parameter record of origin vertex
 left: Face ← NIL,  -- parameter record of face that is to the left of e
 no: NAT      -- edge serial number (for debugging & traversal)
 ];
The routines that manipulate diagrams worry only about maintaining the correction of the edge streucture, as given by the ONext and Dual functions (and their derivatives, like LNext, Sym, Prim, etc.).
-- An interesting special case is when b is dummy, and a is not: the effect will be to add a new loop e interior to the face Right[a] and oriented clockwise. The symmetrical case (a dummy, b normal) is analogous. If both a and b are dummy, the result will be a diagram with a single edge e, and either two faces and one vertex (if a#b) or one face and two vertices (if a=b).
-- The org links of the affected edges will not be fixed by this operator. Therefore, when a vertex is split, the edges coming out of each part will point to the same vertex record; when two vertices are merged into one, the edges coming out of the latter will still point to two different vertex records. The client should use SetOrg[e, v] or SetDest[e, v], as appropriate, after using this operator.
AddNewVertex: PROC [v: Sp.Point, b: Diag.DEdge] RETURNS [e: Diag.DEdge] =
BEGIN
-- Repaints the diagram
 -- To be called with access rights to dr = write
ENABLE UNWIND => {Draw.Release[dr, write]};
Draw.Acquire[dr, write, none];
e ← Diag.ConnectVertices [Diag.MakeSphere[v: NEW [Sp.Point ← v], f: NIL], b];
Draw.DoPaintAll[dr];
Draw.Release[dr, write];
Process.Pause[delay]
END;
AddNewEdge: PROC [a, b: Diag.DEdge] RETURNS [e: Diag.DEdge] =
BEGIN
-- Adds a new edge from Org[a] to Org[b] and inside Left[b]
 -- Also repaints the cube
 -- To be called with access rights to dr = none
ENABLE UNWIND => {Draw.Release[dr, write]};
Draw.Acquire[dr, write, none];
e ← Diag.ConnectVertices [a, b];
Draw.DoPaintAll[dr];
Draw.Release[dr, write];
Process.Pause[delay]
END;
e1 ← cube^ ← AddNewVertex [ [1, 0, 0], cube^ ];
e2 ← AddNewVertex [ [1, 1, 0], e1 ];
e3 ← AddNewVertex [ [0, 1, 0], e2 ];
e4 ← AddNewEdge [ Diag.Sym[e1], e3 ];
e5 ← AddNewVertex [ [0, 1, 1], e3 ];
e6 ← AddNewVertex [ [1, 1, 1], e5 ];
e7 ← AddNewVertex [ [1, 0, 1], e6 ];
e8 ← AddNewVertex [ [0, 0, 1], e7 ];
e9 ← AddNewEdge [ Diag.Sym[e6], e8 ];
e10 ← AddNewEdge [ e2, e6 ];
e11 ← AddNewEdge [ e1, e7 ];
e12 ← AddNewEdge [ e4, e8 ];

(INSERT IN VORONOI - SECOND VERSION - DAMN BUGGED)

erConv: BOOL;
erConv ← FALSE; es ← en;
DO
IF NOT Homo.Convex [p.pt, q.pt, Dest[es].pt] THEN
{IF erConv THEN EXIT}
ELSE
{erConv ← TRUE};
er ← es; es ← Diag.ONext[es]
ENDLOOP;
-- add p to the diagram, connected to q
 ep ← Diag.AddVertex [p, Diag.Sym[es], left];
Bug.out.PutF["\nVorIncr: connected to vertex %g.", int[q.no]];
-- main loop
 DO
-- at this point ep is a new edge out of q, between er and es c.c.w
 r ← Dest[er];
-- delete edges out of q that interfere with p
 DO
et ← Diag.OPrev[er]; t ← Dest[et];
IF NOT Homo.Convex [p.pt, q.pt, t.pt] OR
NOT Homo.Convex [t.pt, r.pt, p.pt] OR
Homo.ShouldConnect24 [p.pt, q.pt, t.pt, r.pt] THEN EXIT;
IF trace THEN ShowEdge [dr: dr, e: er, width: 2, color: blue];
Bug.out.PutF["\nVorIncr: Deleting edge "]; Diag.PutE[er];
[] ← Diag.DeleteEdge [er];
Bug.out.PutF[" - new face: "]; Diag.PutLRing[et];
er ← et; r ← t
ENDLOOP;
-- advance q to next vertex to be connected to p
 DO
es ← Diag.Sym[er]; s ← q;
q ← r; er ← Diag.OPrev[es];
IF q = n OR
NOT Homo.Convex [r.pt, p.pt, s.pt] OR
NOT Homo.Convex [s.pt, q.pt, r.pt] OR
Homo.ShouldConnect24 [r.pt, p.pt, s.pt, q.pt] THEN EXIT;
ENDLOOP;
-- connect p to q (unless q is the starting vertex n) and set Voronoi vertex
 vcur ← MakeVVertex [VVertexCoords[p.pt, s.pt, q.pt]];
IF q # n THEN
{ep ← Diag.Sym[Diag.CloseFace [ep, es, vcur, right]]}
ELSE
{ep ← er; er ← Diag.OPrev[er]; Diag.SetLeft [ep, vcur]};
Bug.out.PutF["\nVorIncr: closed face = "]; Diag.PutLRing [ep];
IF trace THEN
{et ← Diag.LNext [ep];
ShowEdge [dr: dr, e: ep];
DO
et ← Diag.LNext [et];
IF et = ep THEN EXIT;
ShowEdge [dr: dr, e: et, color: red, width: 2]
ENDLOOP};
IF q = n THEN EXIT;
ENDLOOP
END;
VertexPainter, RegionPainter, EdgePainter: Draw.PainterProc;
MakeVertexObject: PROC
[v: REF VertexRec, color: Draw.Color ← Col.black, side: Draw.Size ← 0]
RETURNS [obj: Draw.Object] = INLINE
-- Makes an isolated DVertex or VVertex d into a Drawing object that will paint itself as a square with given side (in pixels).
{RETURN [Draw.Make [VertexPainter, [data: v, color: color, size: side]]]};
MakeRegionObject: PROC
[e: Diag.DEdge, color: Draw.Color ← Col.red, side: Draw.Size ← 0]
RETURNS [obj: Draw.Object] = INLINE
-- The edge e must have a DVertex as origin. The resulting object will normally paint Org[e] as a small square; however, if $VoronoiStyle is in the drawData list, the corresponding Voronoi region will be painted instead.
{RETURN [Draw.Make [RegionPainter, [data: e.rec, color: color, size: side, style: e.ix]]]};
MakeEdgeObject: PROC
[e: Diag.DEdge, color: Draw.Color, width: Draw.Size ← 0]
RETURNS [obj: Draw.Object] = INLINE
-- Makes an isolated Delaunay edge e into a Drawing object that will paint itself in dual or primal form, depending on whether $VoronoiStyle is in the drawData list. The width is in pixels (0 or 1 gives the "thin" line of CedarGraphics).
-- The object's parms.data field points to the whole edge record; parms.style tells which of the four quarters of e is to be considered. It should always be a Delaunay edge (Org[e] :: DVertex), even when $VoronoiStyle is on.
MakeDodecahedron: PROC [dr: Draw.Drawing, dg: REF Diag.Diagram] =
BEGIN OPEN Diag, RealFns;
-- old version - coordinates guessed, not computed
e1, e2, e3, ei1, ei2, ei3, b1, b2, b3, ep1, ep2, ep3: Diag.DEdge;
v0, v1, v2, v3: REF Sp.Point;
f1, f3, ftop, fbot: REF ← NIL;
c, s: REAL;
twopi: REAL = 2.0*3.1415926;
mr: REAL = 1.2; -- guess
 mz: REAL = 0.2; -- guess too
 Draw.Acquire[dr, write, none];
BEGIN ENABLE UNWIND => {Draw.Release[dr, write]};
FOR i: NAT IN [0..5] DO
IF i # 5 THEN
{c ← Cos[twopi*i/5.0];
s ← Sin[twopi*i/5.0];
v0 ← NEW [Sp.Point ← [c, s, -1.0]];
v1 ← NEW [Sp.Point ← [mr*c, mr*s, -mz]];
c ← RealFns.Cos[twopi*(i-0.5)/5.0];
s ← RealFns.Sin[twopi*(i-0.5)/5.0];
v2 ← NEW [Sp.Point ← [mr*c, mr*s, +mz]];
v3 ← NEW [Sp.Point ← [c, s, +1.0]];
e1 ← MakeBridge [v0, v1, NIL];
e2 ← AddVertex [v2, e1, left];
e3 ← AddVertex [v3, e2, left]}
ELSE
{e1 ← ei1; e2 ← ei2; e3 ← ei3};
IF i = 0 THEN
{ei1 ← e1; ei2 ← e2; ei3 ← e3; dg^ ← e1}
ELSE
{b3 ← ConnectVertices [Sym[ep3], Sym[e3]];
b2 ← CloseFace [Sym[e3], ep2, f3, left];
b1 ← CloseFace [Sym[e1], ep1, f1, left];
IF i = 5 THEN
{SetRight [b3, ftop];
SetRight [b1, fbot]}};
ep1 ← e1; ep2 ← e2; ep3 ← e3;
Draw.DoPaintAll[dr];
Process.Pause[delay]
ENDLOOP;
END;
Draw.Release [dr, write]
END;
DrawSegment: PUBLIC PROC
[context: Graphics.Context, sf: Cart.ScaleFactors,
org, dest: Homo.Point, width: Size ← 1, color: Color ← black] = TRUSTED
BEGIN
-- note: size = 0 means standard line size (as defined by DrawTo)
 wd: REAL = Float [width]/2.0;
so, sd: Cart.Point;
[so, sd] ← Homo.ScaleSeg [org, dest, sf];
Graphics.SetColor[context, color];
IF width # 0 THEN
{d: Cart.Vector ← Cart.Sub [sd, so];
tx,ty: REAL;
ds: REAL ← Cart.Length[d];
IF ds > 0.1 THEN {tx ← -wd*d.y/ds; ty ← wd*d.x/ds};
Graphics.MoveTo[context, so.x+tx, so.y+ty];
Graphics.LineTo[context, sd.x+tx, sd.y+ty];
Graphics.LineTo[context, sd.x-tx, sd.y-ty];
Graphics.LineTo[context, so.x-tx, so.y-ty];
Graphics.DrawArea[context]}
ELSE
{Graphics.MoveTo[context, so.x, so.y];
Graphics.LineTo[context, sd.x, sd.y]}
END;
END
ParallelThrough: PROC [r: Line, p: Point] RETURNS [Line];
-- Computes the line parallel to r passing through p. If r is at infinity or indeterminate, or p is indeterminate, or p is at infinity and on the line r, the result is indeterminate. If p is at infinity but not on r, the result is the line at infinity.
ParallelThrough: PROC [r: Line, p: Point] RETURNS [Line];
-- 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.
PerpThrough: PROC [r: Line, p: Point] RETURNS [Line] = INLINE
-- Computes line passing through p and oriented 90 degrees counterclockwise (as seen by p) to the line r. 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 [[r.y*p.w, -r.x*p.w, -r.y*p.x+r.x*p.y]]};
PerpendicularThrough: PROC [v: Cart.Vector, p: Point] RETURNS [Line] = INLINE
-- Computes line passing through p and oriented 90 degrees counterclockwise 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]]};
ParallelThrough: PUBLIC PROC [r: Line, p: Point] RETURNS [Line] =
-- Should be the same as LineThrough [p, Intercept [r, lineAtInfinity]]
{RETURN [r.x*p.w, r.y*p.w, -r.x*p.x-r.y*p.y]};
PerpThrough: PUBLIC PROC [r: Line, p: Point] RETURNS [Line] =
{RETURN [LineThrough [p, ClosestPoint [p, r]]]]};
-- Should be the same as [r.y*p.w, -r.x*p.w, -r.y*p.x+r.x*p.y]
LinePerp: PROC [r: Line] RETURNS [Vector];
-- If the line is at infinity, 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).
LineDir: PROC [r: Line] RETURNS [Cart.Vector];
-- If the line is at infinity, returns the indeterminate vector. Otherwise returns the orientation of r (a unit vector parallel to r as seen from the top plane).
LineNormal: PUBLIC PROC [r: Line] RETURNS [Cart.Vector] =
{mod: REAL = SqRt[r.x*r.x+r.y*r.y];
RETURN [[r.x/mod, r.y/mod]]};
LineDir: PUBLIC PROC [r: Line] RETURNS [Cart.Vector] =
{mod: REAL = SqRt[r.x*r.x+r.y*r.y];
RETURN [[r.y/mod, -r.x/mod]]};
END.