G3dPlaneImpl.mesa
Copyright Ó 1984, 1992 by Xerox Corporation. All rights reserved.
Bloomenthal, October 9, 1992 10:52 am PDT
Andrew Glassner February 19, 1991 4:06 pm PST
Ken Fishkin, November 6, 1991 1:32 pm PST
DIRECTORY G3dBasic, G3dMatrix, G3dPlane, G3dPolygon, G3dVector, Real, RealFns, Rope;
G3dPlaneImpl: CEDAR PROGRAM
IMPORTS G3dMatrix, G3dPolygon, G3dVector, Real, RealFns
EXPORTS G3dPlane
~ BEGIN
Error:     PUBLIC ERROR [reason: ROPE] ~ CODE;
Circle:    TYPE ~ G3dBasic.Circle;
NatSequence:   TYPE ~ G3dBasic.NatSequence;
Pair:     TYPE ~ G3dBasic.Pair;
PairSequence:  TYPE ~ G3dBasic.PairSequence;
PairSequenceRep: TYPE ~ G3dBasic.PairSequenceRep;
Quad:     TYPE ~ G3dBasic.Quad;
Sign:     TYPE ~ G3dBasic.Sign;
Triple:    TYPE ~ G3dBasic.Triple;
TripleSequence:  TYPE ~ G3dBasic.TripleSequence;
Matrix:    TYPE ~ G3dMatrix.Matrix;
MajorPlane:   TYPE ~ G3dPlane.MajorPlane;
Plane:     TYPE ~ G3dPlane.Plane;
PlaneSequence:  TYPE ~ G3dPlane.PlaneSequence;
PlaneSequenceRep: TYPE ~ G3dPlane.PlaneSequenceRep;
Ray:     TYPE ~ G3dPlane.Ray;
ROPE:     TYPE ~ Rope.ROPE;
Plane Creations
FromPolygon: PUBLIC PROC [
points: TripleSequence,
polygon: NatSequence ¬ NIL]
RETURNS [p: Plane]
~ {
normal: Triple ~ G3dPolygon.PolygonNormal[points, polygon, TRUE];
p ¬ [normal.x, normal.y, normal.z, DFromNormalAndPoints[normal, points, polygon], TRUE];
};
FromThreePoints: PUBLIC PROC [p1, p2, p3: Triple, unitize: BOOL ¬ FALSE]
RETURNS [plane: Plane]
~ {
normal: Triple ¬ G3dVector.Cross[G3dVector.Sub[p2, p1], G3dVector.Sub[p3, p1]];
IF G3dVector.Null[normal] THEN ERROR Error["points colinear"];
IF G3dVector.Null[normal] THEN normal ¬ G3dVector.Ortho[G3dVector.Sub[p2, p1]];
plane ¬ FromPointAndNormal[p1, normal, unitize];
plane ¬ [normal.x, normal.y, normal.z,
-(p1.x*(p2.y*p3.z-p3.y*p2.z)-p1.y*(p2.x*p3.z-p3.x*p2.z)+p1.z*(p2.x*p3.y-p3.x*p2.y))];
IF unitize THEN plane ¬ Unit[plane];
};
FromThreePointsAndBehind: PUBLIC PROC [
p1, p2, p3, behind: Triple,
unitize: BOOL ¬ FALSE]
RETURNS [plane: Plane]
~ {
plane ¬ FromThreePoints[p1,p2, p3, unitize];
IF DistanceToPoint[behind, plane] > 0.0 THEN plane ¬ Negate[plane];
};
FromPointAndNormal: PUBLIC PROC [point, normal: Triple, unitize: BOOL ¬ FALSE]
RETURNS [plane: Plane]
~ {
plane ¬ [normal.x, normal.y, normal.z, -normal.x*point.x-normal.y*point.y-normal.z*point.z];
IF unitize THEN plane ¬ Unit[plane];
};
DFromNormalAndPoints: PUBLIC PROC [
normal: Triple,
points: TripleSequence,
polygon: NatSequence ¬ NIL]
RETURNS [REAL]
~ {
nPoints: NAT ¬ IF polygon = NIL THEN points.length ELSE polygon.length;
GetPointN: PROC [n: NAT] RETURNS [Triple] ~ INLINE {
RETURN[IF polygon = NIL THEN points[n] ELSE points[polygon[n]]];
};
min: REAL ¬ 100000.0;
max: REAL ¬ -100000.0;
FOR n: NAT IN [0..nPoints) DO
distance: REAL ¬ G3dVector.Dot[GetPointN[n], normal];
max ¬ MAX[max, distance];
min ¬ MIN[min, distance];
ENDLOOP;
RETURN[-0.5*(min+max)];
};
NormalFromPlane: PUBLIC PROC [plane: Plane, unit: BOOL ¬ FALSE] RETURNS [Triple] ~ {
t: Triple ¬ [plane.x, plane.y, plane.z];
RETURN[IF unit AND NOT plane.normalized THEN G3dVector.Unit[t] ELSE t];
};
Plane Modifications
Unit: PUBLIC PROC [plane: Plane] RETURNS [unitized: Plane] ~ {
IF plane.normalized
THEN RETURN[plane]
ELSE {
length: REAL ¬ G3dVector.Length[[plane.x, plane.y, plane.z]];
IF length = 0.0 THEN ERROR Error["null plane normal"];
IF length = 1.0
THEN unitized ¬ [plane.x, plane.y, plane.z, plane.w, TRUE]
ELSE {
inverse: REAL ¬ 1.0/length;
unitized ¬ [inverse*plane.x, inverse*plane.y, inverse*plane.z, inverse*plane.w];
};
unitized.normalized ¬ TRUE;
};
};
Intersections
IntersectionOf3: PUBLIC PROC [p1, p2, p3: Plane] RETURNS [Triple] ~ {
Intersection is the determinant of a matrix formed from the three planes.
w: REAL ¬ p1.y*(p2.x*p3.z-p3.x*p2.z)-p1.x*(p2.y*p3.z-p3.y*p2.z)-p1.z*(p2.x*p3.y-p3.x*p2.y);
IF w = 0.0 THEN ERROR Error["no intersection"];
RETURN[[
(p1.y*(p2.z*p3.w-p3.z*p2.w)-p1.z*(p2.y*p3.w-p3.y*p2.w)+p1.w*(p2.y*p3.z-p3.y*p2.z))/w,
-(p1.x*(p2.z*p3.w-p3.z*p2.w)-p1.z*(p2.x*p3.w-p3.x*p2.w)+p1.w*(p2.x*p3.z-p3.x*p2.z))/w,
(p1.x*(p2.y*p3.w-p3.y*p2.w)-p1.y*(p2.x*p3.w-p3.x*p2.w)+p1.w*(p2.x*p3.y-p3.x*p2.y))/w]];
};
IntersectionOf2: PUBLIC PROC [plane1, plane2: Plane] RETURNS [line: Ray] ~ {
line.axis ¬ G3dVector.Cross[[plane1.x, plane1.y, plane1.z], [plane2.x, plane2.y, plane2.z]];
IF G3dVector.Null[line.axis] THEN ERROR Error["no intersection"];
line.base ¬ IntersectionOf3[plane1, plane2, [line.axis.x, line.axis.y, line.axis.z, 0.0]];
};
IntersectWithLine: PUBLIC PROC [plane: Plane, line: Ray] RETURNS [Triple] ~ {
alpha: REAL;
pdDot: REAL ¬ G3dVector.Dot[[plane.x, plane.y, plane.z], line.axis];
IF pdDot = 0.0 THEN ERROR Error["no intersection"];
alpha ¬ (-plane.w - G3dVector.Dot[[plane.x, plane.y, plane.z], line.base])/pdDot;
RETURN[G3dVector.ScaleRay[line, alpha]];
};
Projections
ProjectPointToPlane: PUBLIC PROC [point: Triple, plane: Plane] RETURNS [Triple] ~ {
Assumes the plane normal is unit length.
normal: Triple ~ [plane.x, plane.y, plane.z];
distance: REAL ~ G3dVector.Dot[point, normal]+plane.w;
RETURN[G3dVector.Sub[point, G3dVector.Mul[normal, distance]]];
};
ProjectVectorToPlane: PUBLIC PROC [v: Triple, plane: Plane] RETURNS [Triple] ~ {
vv: Triple ¬ G3dVector.Project[v, [plane.x, plane.y, plane.z]];
RETURN[[v.x-vv.x, v.y-vv.y, v.z-vv.z]];
};
ProjectPointToMajorPlane: PUBLIC PROC [p: Triple, majorPlane: MajorPlane]
RETURNS [Pair]
~ {
Note: Ned Greene claims to should be called "principal" planes.
Jules Bloomenthal notes:
Follow me for a sec: Algorithm to compute intersection of line with triangle:
 compute plane of triangle
 compute major plane of triangle by finding which component of plane normal is
greatest
  if z is greatest, major plane is xy, if y, then xz, if x, then yz
 project triangle vertices to major plane
 compute 2d line equations of these projected points

 now, compute intersection of line with plane and project intersection to
 major plane. now determine if projection is inside triangle by
 multiplying 2d point by 2d line; if all results are positive, pt is in
triangle

Or so I thought. There are two subtleties:
 1. the computation of the 2d line equation depends on the ordering of the
 triangle vertices. This ordering can be determined from the sign of the
 component of greatest magnitude. ok, easy fix.

 2. a projection to the major plane has an implied ordering. we're going from
 (x,y,z) to (x,y). If the major plane is xy, then (x,y,z) -> (x,y), if the
major plane
 is xz, then (x,y,z) -> (x,z); and if the major plane is yz, then (x,y,z) ->
(y,z).
 Right?
 Wrong!
 Because the axes are implicitly ordered x, y, z, if the major plane is xz,
then
 (x,y,z) should --> (z, x)!! Right?? I think so. Do you agree?
 Anyway, this is a subtlety easily overlooked. I just wasted about two days
 trying to figure out what was going wrong in a complex 3d algorithm that
 depended at a fairly low level on the correct intersection of a line with
 a triangle.

Question: might this subtlety have a name? Something like the "asymmetric
nature
of 3d to 2d projection?" Or is the asymmetry really an illusion that comes
because
we blindly project (x,y,z) to (x,z) when it should be (z, x)?
Paul Heckbert replies:
Yes, this is an interesting subtlety.
Related is the formula for the cross product of two vectors.
Those formulas are:
 (AxB)1 = a2*b3 - a3*b2
 (AxB)2 = a3*b1 - a1*b3
 (AxB)3 = a1*b2 - a2*b1

where A=(a1,a2,a3) and B=(b1,b2,b3).
Note that the second one is a3b1-a1b3 and not a1b3-a3b1 as you might
think if you sort the terms (the string "a1b3" comes before "a3b1"
lexicographically).
The way I remember this formula is to remember the first term of the
first line "(AxB)1 = a2*b3 - ..." and then remember that this is a
2x2 determinant (a "minor" of a 3x3 matrix) and 2x2 determinants
always involve products of the form foo[i]*bar[j]-foo[j]*bar[i].
So that gives me the first line.

 (AxB)1 = a2*b3 - a3*b2
 (AxB)? = a?*b? - a?*b?
 (AxB)? = a?*b? - a?*b?

Then I cyclically permute the numbers as I go down each column of unknown '?'s,
keeping the letters the same from line to line.
If you start at 1, then the cycle is 1,2,3.
Starting at 2, the cycle is 2,3,1; and if you start at 3, the cycle is 3,1,2:

    \
   1 -------@ 2
    @-- / /
    |\ /
    \ |/
    \ @--
    3

Similarly, if you've got orthonormal (orthogonal and multually perpendicular)
basis vectors I, J, and K where K=IxJ, (e.g. I=(1,0,0), J=(0,1,0), K=(0,0,1)),
then the rule is that I=JxK and J=KxI. If you line them up, the cycles
are easier to see:

 K = I x J
 I = J x K
 J = K x I

As you go down each column, this is a clockwise cycle of the symbols (I,J,K).
Since I, J, and K are consecutive letters of the alphabet, you can also
remember the rules by using C-like reasoning: increment or decrement the
char's by 1, with wraparound, and you get a self-consistent set of rules.
If you used a counterclockwise cycle, you'd have a self-consistent set of
cross-product rules for a left-handed coordinate frame.

When you're doing this projection of 3-D to 2-D business, it helps
to use the same kind of cyclic rule to avoid throwing ugly minus signs
into the formulas:

(In the cross product formulas, you could always just write them
 (AxB)1 = a2*b3 - a3*b2
 (AxB)2 = -a1*b3 + a3*b1
 (AxB)3 = a1*b2 - a2*b1
but the code doesn't line up so nicely or look so pretty,
so I never write the formulas this way.

As for a name for this phenomenon, I'd call it the cyclic nature of
cross products in 3-D.

-Paul

ps: Haines' chapter in "Introduction to Ray Tracing", Glassner, ed,
has ray-polygon intersection formulas. He gets things right, but
I don't think his derivations are very elegant.
RETURN[SELECT majorPlane FROM
xy =>   [p.x, p.y],
yz =>   [p.y, p.z],
ENDCASE => [p.z, p.x]]; -- was [p.x, p.z]
};
ProjectPointsToXYPlane: PUBLIC PROC [points: TripleSequence]
RETURNS [pairs: PairSequence]
~ {
pairs ¬ NEW[PairSequenceRep[points.length]];
pairs.length ¬ points.length;
FOR n: NAT IN [0..points.length) DO
pairs[n] ¬ [points[n].x, points[n].y];
ENDLOOP;
};
ProjectPointsToXZPlane: PUBLIC PROC [points: TripleSequence]
RETURNS [pairs: PairSequence]
~ {
pairs ¬ NEW[PairSequenceRep[points.length]];
pairs.length ¬ points.length;
FOR n: NAT IN [0..points.length) DO
pairs[n] ¬ [points[n].x, points[n].z];
ENDLOOP;
};
ProjectPointsToYZPlane: PUBLIC PROC [points: TripleSequence]
RETURNS [pairs: PairSequence]
~ {
pairs ¬ NEW[PairSequenceRep[points.length]];
pairs.length ¬ points.length;
FOR n: NAT IN [0..points.length) DO
pairs[n] ¬ [points[n].y, points[n].z];
ENDLOOP;
};
ProjectPointsToMajorPlane: PUBLIC PROC [points: TripleSequence, majorPlane: MajorPlane]
RETURNS [pairs: PairSequence]
~ {
pairs ¬ SELECT majorPlane FROM
xy => ProjectPointsToXYPlane[points],
xz => ProjectPointsToXZPlane[points],
ENDCASE => ProjectPointsToYZPlane[points];
};
Miscellany
Side: PUBLIC PROC [point: Triple, plane: Plane] RETURNS [Sign] ~ {
s: REAL ¬ point.x*plane.x+point.y*plane.y+point.z*plane.z+plane.w;
RETURN[IF s > 0.0 THEN positive ELSE negative];
};
Negate: PUBLIC PROC [plane: Plane] RETURNS [Plane] ~ {
RETURN[[-plane.x, -plane.y, -plane.z, -plane.w]];
};
AlignWithXYPlane: PUBLIC PROC [plane: Plane, out: Matrix ¬ NIL] RETURNS [Matrix] ~ {
v: Triple ¬ [plane.x, plane.y, plane.z];
sqLen: REAL ¬ G3dVector.Square[v];
normal: Triple ¬ IF sqLen = 1.0 THEN v ELSE G3dVector.Div[v, RealFns.SqRt[sqLen]];
axis: Triple ¬ G3dVector.Cross[normal, G3dBasic.zAxis];
degrees: REAL ¬ G3dVector.AngleBetween[normal, G3dBasic.zAxis];
center: Triple ¬ CenterOfPlane[plane];
out ¬ G3dMatrix.MakeRotate[axis, degrees, TRUE, out];
out ¬ G3dMatrix.Translate[out, G3dVector.Negate[center], out];
RETURN[out];
};
ClosestPair: PUBLIC PROC [pairs: PairSequence, pair: Pair] RETURNS [index: NAT] ~ {
min: REAL ¬ 1000000.0;
IF pairs = NIL THEN RETURN[0];
FOR n: NAT IN [0..pairs.length) DO
p: Pair ~ pairs[n];
dx: REAL ~ p.x-pair.x;
dy: REAL ~ p.y-pair.y;
d: REAL ~ dx*dx+dy*dy;
IF d < min THEN {index ¬ n; min ¬ d};
ENDLOOP;
};
GetMajorPlane: PUBLIC PROC [plane: Plane] RETURNS [MajorPlane] ~ {
x: REAL ~ ABS[plane.x];
y: REAL ~ ABS[plane.y];
z: REAL ~ ABS[plane.z];
RETURN[SELECT MAX[x, y, z] FROM
x   => yz,
y   => xz,
ENDCASE => xy];
};
CenterOfPlane: PUBLIC PROC [plane: Plane] RETURNS [Triple] ~ {
s: REAL ~ -plane.w/(plane.x*plane.x+plane.y*plane.y+plane.z*plane.z);
RETURN[[s*plane.x, s*plane.y, s*plane.z]];
};
DistanceToPoint: PUBLIC PROC [point: Triple, plane: Plane, unitize: BOOL ¬ TRUE]
RETURNS [d: REAL]
~ {
d ¬ point.x*plane.x+point.y*plane.y+point.z*plane.z+plane.w;
IF unitize THEN {
length: REAL ¬ G3dVector.Length[[plane.x, plane.y, plane.z]];
IF length = 0.0 THEN ERROR Error["null plane normal"];
d ¬ d/length;
};
};
CircleFromPoints: PUBLIC PROC [p1, p2, p3: Triple] RETURNS [circle: Circle] ~ {
plane: Plane ¬ FromThreePoints[p1, p2, p3];
v0: Triple ¬ G3dVector.Unit[
G3dVector.Cross[G3dVector.Sub[p1, p2], [plane.x, plane.y, plane.z]]];
v1: Triple ¬ G3dVector.Unit[
G3dVector.Cross[G3dVector.Sub[p2, p3], [plane.x, plane.y, plane.z]]];
dot: REAL ¬ G3dVector.Dot[v0, v1];
IF dot = 1.0
THEN ERROR Error["points colinear"]
ELSE {
m0: Triple ¬ G3dVector.Midpoint[p1, p2];
m1: Triple ¬ G3dVector.Midpoint[p2, p3];
dm: Triple ¬ G3dVector.Sub[m1, m0];
t: REAL ¬
(dot*G3dVector.Dot[v0, dm]-G3dVector.Dot[v1, dm])/(1.0-dot*dot);
circle.center ¬ G3dVector.ScaleRay[[m1, v1], t];
circle.radius ¬ G3dVector.Distance[circle.center, p1];
};
};
Plane Sequences
CopyPlaneSequence: PUBLIC PROC [planes: PlaneSequence] RETURNS [PlaneSequence] ~ {
copy: PlaneSequence ¬ NIL;
IF planes # NIL THEN {
copy ¬ NEW[PlaneSequenceRep[planes.length]];
copy.length ¬ planes.length;
FOR n: NAT IN [0..planes.length) DO copy[n] ¬ planes[n]; ENDLOOP;
};
RETURN[copy];
};
AddToPlaneSequence: PUBLIC PROC [planes: PlaneSequence, plane: Plane] RETURNS [PlaneSequence]
~ {
IF planes = NIL THEN planes ¬ NEW[PlaneSequenceRep[1]];
IF planes.length = planes.maxLength THEN planes ¬ LengthenPlaneSequence[planes];
planes[planes.length] ¬ plane;
planes.length ¬ planes.length+1;
RETURN[planes];
};
LengthenPlaneSequence: PUBLIC PROC [planes: PlaneSequence, amount: REAL ¬ 1.3]
RETURNS [new: PlaneSequence] ~ {
newLength: NAT ¬ MAX[Real.Ceiling[amount*planes.maxLength], 3];
new ¬ NEW[PlaneSequenceRep[newLength]];
FOR i: NAT IN [0..planes.length) DO new[i] ¬ planes[i]; ENDLOOP;
new.length ¬ planes.length;
};
END.