-- File: SVPolygon3dImpl.mesa
-- Last edited by Bier on December 18, 1982 1:11 am
-- Author: Eric Bier on August 21, 1982 1:49 pm
-- Contents: Most of the code necessary to describe the behavior of planes and polygons in space. (including the procedures necessary to make planar polygons act as master objects). Part of the Solidviews 3d Illustrator.
DIRECTORY
Matrix3d,
RealFns,
SVPolygon2d,
SVPolygon3d,
SVVector3d;
SVPolygon3dImpl: PROGRAM
IMPORTS RealFns, SVPolygon2d
EXPORTS SVPolygon3d =
BEGIN
Plane: TYPE = REF PlaneObj;
PlaneObj: TYPE = SVPolygon3d.PlaneObj;
Point3d: TYPE = Matrix3d.Point3d;
Poly3d: TYPE = REF Poly3dObj;
Poly3dObj: TYPE = SVPolygon3d.Poly3dObj;
Polygon: TYPE = REF PolygonObj;
PolygonObj: TYPE = SVPolygon2d.PolygonObj;
Vector: TYPE = SVVector3d.Vector;
-- GLOBAL VARIABLES
root3, root3Over3, twoRoot3Over3: REAL; -- constants computed at load time in Init (see below).
CreatePoly: PUBLIC PROC [len: NAT] RETURNS [poly3d: Poly3d] = {
poly3d ← NEW[Poly3dObj[len]];
poly3d.len ← 0;
};
CircumHexagon: PUBLIC PROC [y: REAL, r: REAL] RETURNS [hex: Poly3d] = {
-- creates a hexigon which circumscribes the circle (in the Y = y plane) whose center is at [0, y, 0] and whose radius is r.
hex ← CreatePoly[6];
hex ← AddPolyPoint[hex, [r*root3Over3, y, r]]; -- front right
hex ← AddPolyPoint[hex, [r*twoRoot3Over3, y, 0]]; -- right
hex ← AddPolyPoint[hex, [r*root3Over3, y, -r]]; -- back right
hex ← AddPolyPoint[hex, [-r*root3Over3, y, -r]]; -- back left
hex ← AddPolyPoint[hex, [-r*twoRoot3Over3, y, 0]]; -- left
hex ← AddPolyPoint[hex, [-r*root3Over3, y, r]]; -- front left
};
ClearPoly: PUBLIC PROC [poly: Poly3d] = {
poly.len ← 0;
};
AddPolyPoint: PUBLIC PROC [poly: Poly3d, point: Point3d] RETURNS [polyPlusPoint: Poly3d] = {
IF poly.len = poly.maxVerts THEN {
polyPlusPoint ← CreatePoly[poly.maxVerts+1];
polyPlusPoint ← PartPolyGetsPartPoly[poly, 0, polyPlusPoint, 0, poly.maxVerts];
polyPlusPoint[poly.maxVerts] ← point;
polyPlusPoint.len ← poly.maxVerts + 1;
}
ELSE {
poly[poly.len] ← point;
poly.len ← poly.len + 1;
polyPlusPoint ← poly;
};
}; -- end of AddPolyPoint
PutPolyPoint: PUBLIC PROC [poly: Poly3d, index: NAT, point: Point3d] RETURNS [newPoly: Poly3d] = {
IF index+1 > poly.maxVerts THEN
{newPoly ← CreatePoly[index+1];
newPoly ← PartPolyGetsPartPoly[poly, 0, newPoly, 0, poly.maxVerts];
newPoly[index] ← point;
newPoly.len ← index+1;
}
ELSE {
IF index+1 > poly.len THEN poly.len ← index+1;
poly[index] ← point;
newPoly ← poly};
};
PartPolyGetsPartPoly: PUBLIC PROC [fromPoly: Poly3d, fromStart: NAT, toPoly: Poly3d, toStart: NAT, duration: NAT] RETURNS [newPoly: Poly3d] = {
-- toPoly in range [toStart..toStart+duration-1] gets fromPoly in range [fromStart..fromStart+duration-1];
j: NAT ← fromStart;
IF toStart+duration > toPoly.maxVerts THEN {
newPoly ← CreatePoly[toStart+duration];
newPoly ← PartPolyGetsPartPoly[toPoly, 0, newPoly, 0, toPoly.len];
newPoly.len ← toStart+duration}
ELSE {newPoly ← toPoly;
IF toStart+duration > newPoly.len THEN newPoly.len ← toStart+duration};
IF fromStart+duration > fromPoly.len THEN ERROR; -- don't use bad data
FOR i: NAT IN[toStart..toStart+duration-1] DO
newPoly[i] ← fromPoly[j];
j ← j + 1;
ENDLOOP;
};
PlaneFromPoly3d: PUBLIC PROC [poly: Poly3d] RETURNS [plane: Plane] = {
-- assumes that the normal of the plane should point in the direction such if the normal points toward a viewer, the viewer will perceive the points of poly to run clockwise.
-- This algorithm was taken from Foley and Van Dam "Fundamentals of Interactive Computer Graphics" p. 512-513, where the authors note that for the plane equation Ax + By + Cz + D = 0,
-- A, B, and C have the same relationship to each other as the areas of the polygon when it is projected onto the yz, xz, and xy planes respectively. Hence we can just set A, B, and C equal to these respective values and then plug in a point to solve for D. The areas are computed by the signed areas of the trapezoids formed by each each and its projection onto the y, z, and x axis respectively for A, B, and C
iPlusOne: NAT;
iPoint, iPlusOnePoint: Point3d;
plane ← NEW[PlaneObj];
plane.A ← 0;
plane.B ← 0;
plane.C ← 0;
FOR i: NAT IN[1..poly.len] DO
iPlusOne ← IF i = poly.len THEN 1 ELSE i + 1;
iPoint ← poly[i - 1];
iPlusOnePoint ← poly[iPlusOne - 1];
plane.A ← plane.A + (iPoint[3] + iPlusOnePoint[3])*(iPlusOnePoint[2] - iPoint[2]);
plane.B ← plane.B + (iPoint[1] + iPlusOnePoint[1])*(iPlusOnePoint[3] - iPoint[3]);
plane.C ← plane.C + (iPoint[2] + iPlusOnePoint[2])*(iPlusOnePoint[1] - iPoint[1]);
ENDLOOP;
-- now find D. Use poly[0] as the point to substitute in.
iPoint ← poly[0];
plane.D ← -(iPoint[1]*plane.A + iPoint[2]*plane.B + iPoint[3]*plane.C);
}; -- end of PlaneFromPoly3d
PlaneFromPointAndNormal: PUBLIC PROC [point: Point3d, normal: Vector] RETURNS [plane: Plane] = {
plane ← NEW[PlaneObj];
plane.A ← normal[1];
plane.B ← normal[2];
plane.C ← normal[3];
-- Ax + By + Cz + D = 0. D = -(Ax + By + Cz)
plane.D ← -(plane.A*point[1] + plane.B*point[2] + plane.C*point[3]);
}; -- end of PlaneFromPointAndNormal
NormalOfPlane: PUBLIC PROC [plane: Plane] RETURNS [normal: Vector] = {
normal[1] ← plane.A;
normal[2] ← plane.B;
normal[3] ← plane.C;
}; -- end of NormalFromPlane
MaximumDeviationWithinPoly3d: PROC [poly: Poly3d] RETURNS [distance: REAL] = {
-- finds the distance of each point in the poly from the plane of the poly and returns the maximum such distance
plane: Plane ← PlaneFromPoly3d[poly];
SIGNAL NotYetImplemented; -- **** not yet implemented
};
NotYetImplemented: SIGNAL = CODE;
SignedPointToPlaneDistance: PUBLIC PROC [point: Point3d, plane: Plane] RETURNS [distance: REAL] = {
distance ← (plane.A*point[1] + plane.B*point[2] + plane.C*point[3] + plane.D)/
RealFns.SqRt[plane.A*plane.A + plane.B*plane.B + plane.C*plane.C];
}; -- end of SignedPointToPlaneDistance
PointOnNormalSideOfPlane: PUBLIC PROC [point: Point3d, plane: Plane] RETURNS [BOOL] = {
plugInPoint: REAL ← (plane.A*point[1] + plane.B*point[2] + plane.C*point[3] + plane.D);
RETURN [plugInPoint > 0];
};
Quadrant: TYPE = SVPolygon3d.Quadrant; -- {inin, inout, outin, outout};
QuadrantOfPoint: PUBLIC PROC [point: Point3d, plane1, plane2: Plane] RETURNS [Quadrant] = {
-- classifys point with respect to two planes. outin means the point is on the normal side of the first plane and the non-normal side of the other plane, etc.
IF PointOnNormalSideOfPlane[point, plane1] THEN {
IF PointOnNormalSideOfPlane[point, plane2] THEN RETURN[outout]
ELSE RETURN[outin]}
ELSE {
IF PointOnNormalSideOfPlane[point, plane2] THEN RETURN[inout]
ELSE RETURN[inin]};
}; -- QuadrantOfPoint
ProjectPolyToXYPlane: PUBLIC PROC [poly3d: Poly3d] RETURNS [polygon: Polygon] = {
polygon ← SVPolygon2d.CreatePoly[poly3d.len];
FOR i: NAT IN[0..poly3d.len) DO
polygon[i][1] ← poly3d[i][1];
polygon[i][2] ← poly3d[i][2];
ENDLOOP;
polygon.len ← poly3d.len;
};
ProjectPolyToYZPlane: PUBLIC PROC [poly3d: Poly3d] RETURNS [polygon: Polygon] = {
FOR i: NAT IN[0..poly3d.len) DO
polygon[i][1] ← poly3d[i][2];
polygon[i][2] ← poly3d[i][3];
ENDLOOP;
polygon.len ← poly3d.len;
};
ProjectPolyToXZPlane: PUBLIC PROC [poly3d: Poly3d] RETURNS [polygon: Polygon] = {
FOR i: NAT IN[0..poly3d.len) DO
polygon[i][1] ← poly3d[i][1];
polygon[i][2] ← poly3d[i][3];
ENDLOOP;
polygon.len ← poly3d.len;
};
Init: PROC = {
root3 ← RealFns.SqRt[3];
root3Over3 ← root3/3;
twoRoot3Over3 ← 2*root3Over3;
};
Init[];
END.