<<>> <> <> <> <> <> 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; <> 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 normal ¬ G3dVector.Ortho[G3dVector.Sub[p2, p1]]; plane ¬ FromPointAndNormal[p1, normal, unitize]; <> <<-(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))];>> <> }; 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]; }; <> 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; }; }; <> IntersectionOf3: PUBLIC PROC [p1, p2, p3: Plane] RETURNS [Triple] ~ { <> 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]]; }; <> ProjectPointToPlane: PUBLIC PROC [point: Triple, plane: Plane] RETURNS [Triple] ~ { <> 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] ~ { <> <> < (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]; }; <> 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]; }; }; <> 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.