Polygons3dImpl.mesa
Copyright © 1985 by Xerox Corporation. All rights reserved.
Bloomenthal, February 26, 1987 7:27:20 pm PST
DIRECTORY Polygons3d, Real, Vector2, Vector3d;
Polygons3dImpl: CEDAR PROGRAM
IMPORTS Real, Vector2, Vector3d
EXPORTS Polygons3d
~ BEGIN
OPEN Polygons3d;
Normal Procedures
TriangleNormal: PUBLIC PROC [p0, p1, p2: Triple] RETURNS [Triple] ~ {
RETURN[[Vector3d.Cross[Vector3d.Sub[p1, p0], Vector3d.Sub[p2, p1]]];
v1: Triple ← [p1.x-p0.x, p1.y-p0.y, p1.z-p0.z];
v2: Triple ← [p2.x-p1.x, p2.y-p1.y, p2.z-p1.z];
cross: Triple ← [v1.y*v2.z-v1.z*v2.y, v1.z*v2.x-v1.x*v2.z, v1.x*v2.y-v1.y*v2.x];
m: REAL ← cross.x*cross.x+cross.y*cross.y+cross.z*cross.z;
RETURN[IF m # 0.0 THEN [cross.x/m, cross.y/m, cross.z/m] ELSE cross];
RETURN[[
(p0.y*p1.z + p1.y*p2.z + p2.y*p0.z - p0.y*p2.z - p1.y*p0.z - p2.y*p1.z),
(p0.z*p1.x + p1.z*p2.x + p2.z*p0.x - p0.z*p2.x - p1.z*p0.x - p2.z*p1.x),
(p0.x*p1.y + p1.x*p2.y + p2.x*p0.y - p0.x*p2.y - p1.x*p0.y - p2.x*p1.y)]];
};
PolygonNormal: PUBLIC PROC [vertices: TripleSequence] RETURNS [Triple] ~ {
normal: Triple ← [0.0, 0.0, 0.0];
FOR i: NAT IN[0..vertices.length) DO        -- Newells' method
p0: Triple ← vertices[i];
p1: Triple ← vertices[(i+1) MOD vertices.length];
normal ← [
normal.x+(p0.y-p1.y)*(p0.z+p1.z),
normal.y+(p0.z-p1.z)*(p0.x+p1.x),
normal.z+(p0.x-p1.x)*(p0.y+p1.y)];
ENDLOOP;
RETURN[normal];
};
PolygonNormals: PUBLIC PROC [
polygons: REF NatTable,
vertices: TripleSequence,
normals: TripleSequence ← NIL]
RETURNS [TripleSequence]
~ {
IF polygons # NIL AND vertices # NIL THEN {
IF normals = NIL OR normals.length < polygons.length
THEN normals ← NEW[TripleSequenceRep[polygons.length]];
normals.length ← polygons.length;
FOR n: NAT IN [0..polygons.length) DO
poly: REF NatSequence ← polygons[n];
nSides: INTEGER ← poly.length;
normals[n] ← [0.0, 0.0, 0.0];
FOR i: NAT IN[0..nSides) DO        -- Newells' method
p0: Triple ← vertices[poly[i]];
p1: Triple ← vertices[poly[(i+1) MOD nSides]];
normals[n] ← [
normals[n].x+(p0.y-p1.y)*(p0.z+p1.z),
normals[n].y+(p0.z-p1.z)*(p0.x+p1.x),
normals[n].z+(p0.x-p1.x)*(p0.y+p1.y)];
ENDLOOP;
ENDLOOP;
};
RETURN[normals];
};
ApplyToFrontFacingPolygons: PUBLIC PROC [
proc: PROC[nPoly: NAT],
polygons: REF NatTable,
pairs: PairSequence,
view: Matrix,
normals: TripleSequence]
~ {
IF polygons = NIL OR pairs = NIL OR view = NIL OR normals = NIL THEN RETURN;
FOR n: NAT IN[0..polygons.length) DO
IF normals[n].x*view[0][2]+normals[n].y*view[1][2]+normals[n].z*view[2][2] >= 0.0
THEN proc[n];
ENDLOOP;
};
Miscellaneous Procedures
CopyPolygon: PUBLIC PROC [poly: REF NatSequence] RETURNS [copy: REF NatSequence] ~ {
copy ← NEW[NatSequence[poly.length]];
copy.length ← poly.length;
FOR n: NAT IN [0..poly.length) DO
copy[n] ← poly[n];
ENDLOOP;
};
ReversePolygon: PUBLIC PROC [poly: REF NatSequence] RETURNS [REF NatSequence] ~ {
FOR n: NAT IN [0..poly.length/2) DO
nn: NAT ← poly.length-n-1;
temp: NAT ← poly[n];
poly[n] ← poly[nn];
poly[nn] ← temp;
ENDLOOP;
RETURN[poly];
};
CenterOfPolygon: PUBLIC PROC [poly: REF NatSequence, points: TripleSequence]
RETURNS [v: Triple ← origin] ~ {
FOR n: NAT IN [0..poly.length) DO
v ← Vector3d.Add[v, points[poly[n]]];
ENDLOOP;
IF poly.length # 0 THEN v ← Vector3d.Div[v, poly.length];
};
CentersOfPolygons: PUBLIC PROC [
polygons: REF NatTable,
points: TripleSequence,
centers: TripleSequence ← NIL]
RETURNS [TripleSequence]
~ {
IF polygons = NIL OR points = NIL THEN RETURN[NIL];
IF centers = NIL OR centers.maxLength < polygons.length
THEN centers ← NEW[TripleSequenceRep[polygons.length]];
centers.length ← polygons.length;
FOR n: NAT IN [0..polygons.length) DO
centers[n] ← CenterOfPolygon[polygons[n], points];
ENDLOOP;
RETURN[centers];
};
ProjectToXYPlane: PUBLIC PROC [
poly: REF NatSequence, points: TripleSequence, o, x, y: Triple]
RETURNS [pairs: PairSequence]
~ {
pairs ← NEW[PairSequenceRep[poly.length]];
pairs.length ← poly.length;
FOR n: NAT IN [0..pairs.length) DO
q: Triple ~ Vector3d.Sub[points[poly[n]], o];
pairs[n] ← [Vector3d.Dot[q, x], Vector3d.Dot[q, y]];
ENDLOOP;
};
GetMinMaxPairs: PUBLIC PROC [pairs: PairSequence] RETURNS [mm: MinMaxPairs] ~ {
mm ← [[10000.0, 10000.0], [-10000.0, -10000.0]];
FOR n: NAT IN [0..pairs.length) DO
p: Pair ~ pairs[n];
mm.min ← [MIN[mm.min.x, p.x], MIN[mm.min.y, p.y]];
mm.max ← [MAX[mm.max.x, p.x], MAX[mm.max.y, p.y]];
ENDLOOP;
};
GetMinMaxTriples: PUBLIC PROC [triples: TripleSequence] RETURNS [mm: MinMaxTriples] ~ {
mm ← [[10000.0, 10000.0, 10000.0], [-10000.0, -10000.0, -10000.0]];
FOR n: NAT IN [0..triples.length) DO
p: Triple ~ triples[n];
mm.min ← [MIN[mm.min.x, p.x], MIN[mm.min.y, p.y], MIN[mm.min.z, p.z]];
mm.max ← [MAX[mm.max.x, p.x], MAX[mm.max.y, p.y], MAX[mm.max.z, p.z]];
ENDLOOP;
};
SquarePairDistance: PUBLIC PROC [p0, p1: Pair] RETURNS [REAL] ~ {
RETURN[Vector2.Square[Vector2.Sub[p0, p1]]];
};
SquareTripleDistance: PUBLIC PROC [p0, p1: Triple] RETURNS [REAL] ~ {
RETURN[Vector3d.SquareLength[Vector3d.Sub[p0, p1]]];
};
Triangle Procedures
LengthenTriangleSequence: PUBLIC PROC [t: TriangleSequence]
RETURNS [TriangleSequence] ~ {
old: TriangleSequence ~ t;
t ← NEW[TriangleSequenceRep[Real.RoundI[1.3*old.length]]];
FOR i: NAT IN [0..old.length) DO t[i] ← old[i]; ENDLOOP;
t.length ← old.length;
RETURN[t];
};
Triangulate2Polygons: PUBLIC PROC [
poly0, poly1: REF NatSequence, points: TripleSequence]
RETURNS [triangles: TriangleSequence]
~ {
NewTriangle: PROC [a, b, c: NAT] ~ {
IF triangles.length >= triangles.maxLength
THEN triangles ← LengthenTriangleSequence[triangles];
triangles[triangles.length] ← [a, b, c];
triangles.length ← triangles.length+1;
};
ScaleAndCenter: PROC ~ {
PutInUnitSquare: PROC [pairs: PairSequence] ~ {
FOR n: NAT IN [0..pairs.length) DO
pairs[n] ← [sx*(pairs[n].x-cx), sy*(pairs[n].y-cy)];
ENDLOOP;
};
amm: MinMaxPairs ← GetMinMaxPairs[apairs];
bmm: MinMaxPairs ← GetMinMaxPairs[bpairs];
mm: MinMaxPairs ← [
[MIN[amm.min.x, bmm.min.x], MIN[amm.min.y, bmm.min.y]],
[MAX[amm.max.x, bmm.max.x], MAX[amm.max.y, bmm.max.y]]];
cx: REAL ← 0.5*(mm.min.x+mm.max.x);
cy: REAL ← 0.5*(mm.min.y+mm.max.y);
sx: REAL ← IF mm.min.x # mm.max.x THEN 2.0/(mm.max.x-mm.min.x) ELSE 0.0;
sy: REAL ← IF mm.min.y # mm.max.y THEN 2.0/(mm.max.y-mm.min.y) ELSE 0.0;
PutInUnitSquare[apairs];
PutInUnitSquare[bpairs];
};
SetStartingPoints: PROC ~ {
s, max: REAL ← 10000.0;
FOR n: NAT IN [0..a.length) DO
FOR nn: NAT IN [0..b.length) DO
IF (s ← SquarePairDistance[apairs[n], bpairs[nn]]) < max
THEN {astart ← a0 ← n; bstart ← b0 ← nn; max ← s};
ENDLOOP;
ENDLOOP;
};
a: REF NatSequence ~ poly0;
b: REF NatSequence ~ ReversePolygon[CopyPolygon[poly1]];
acenter: Triple ← CenterOfPolygon[a, points];
bcenter: Triple ← CenterOfPolygon[b, points];
center: Triple ← Vector3d.Mul[Vector3d.Add[acenter, bcenter], 0.5];
z: Triple ← Vector3d.Normalize[Vector3d.Sub[bcenter, acenter]];
x: Triple ← Vector3d.Ortho[z, yAxis];
y: Triple ← Vector3d.Normalize[Vector3d.Cross[x, z]];
apairs: PairSequence ← ProjectToXYPlane[a, points, center, x, y];
bpairs: PairSequence ← ProjectToXYPlane[b, points, center, x, y];
a0, a1, astart, b0, b1, bstart: NAT ← 0;
nTriangles: NAT ← 0;
maxNTriangles: NAT ← poly0.length+poly1.length+2; -- +2 is for good luck
ScaleAndCenter[];
SetStartingPoints[];
a1 ← (a0+1) MOD a.length;
b1 ← (b0+1) MOD b.length;
triangles ← NEW[TriangleSequenceRep[points.length]];
DO
IF SquarePairDistance[apairs[a0], bpairs[b1]] < SquarePairDistance[apairs[a1], bpairs[b0]]
THEN {
NewTriangle[a[a0], b[b0], b[b1]];
b0 ← b1;
b1 ← (b1+1) MOD b.length;
}
ELSE {
NewTriangle[a[a0], b[b0], a[a1]];
a0 ← a1;
a1 ← (a1+1) MOD a.length;
};
IF a0 = astart AND b0 = bstart THEN EXIT;
IF (nTriangles ← nTriangles+1) >= maxNTriangles THEN EXIT;
ENDLOOP;
};
END.