DIRECTORY Algebra3d, CSG, Matrix3d, SVFaces, SVFacesCast, SVVector3d; SVFacesCastImpl: PROGRAM IMPORTS Algebra3d, SVVector3d EXPORTS SVFacesCast = BEGIN Cone: TYPE = REF ConeObj; ConeObj: TYPE = SVFaces.ConeObj; DiskRing: TYPE = REF DiskRingObj; DiskRingObj: TYPE = SVFaces.DiskRingObj; Cylinder: TYPE = REF CylinderObj; CylinderObj: TYPE = SVFaces.CylinderObj; EdgeOnRect: TYPE = REF EdgeOnRectObj; EdgeOnRectObj: TYPE = SVFaces.EdgeOnRectObj; Point3d: TYPE = Matrix3d.Point3d; Ray: TYPE = REF RayObj; RayObj: TYPE = CSG.RayObj; Vector: TYPE = SVVector3d.Vector; FaceHitRec: TYPE = REF FaceHitRecObj; FaceHitRecObj: TYPE = SVFacesCast.FaceHitRecObj; almostZero: REAL = 1.0e-12; Abs: PROC [r: REAL] RETURNS [REAL] = { RETURN [IF r >= 0 THEN r ELSE -r]; }; PositiveRoots: PROCEDURE [a, b, c: REAL] RETURNS [rootArray: ARRAY[1..2] OF REAL, rootCount: NAT] = { i: NAT _ 1; [rootArray, rootCount] _ Algebra3d.QuadraticFormula[a, b, c]; IF rootCount = 0 THEN RETURN; WHILE i <= rootCount DO IF rootArray[i] <= 0 THEN { FOR j: NAT IN [i..rootCount-1] DO rootArray[j] _ rootArray[j+1]; ENDLOOP; rootCount _ rootCount - 1; } ELSE {i _ i + 1}; ENDLOOP; }; ConeCast: PUBLIC PROC [cone: Cone, localRay: Ray] RETURNS [hits: FaceHitRec] = { a, b, c: REAL; y: REAL; roots: ARRAY[1..2] OF REAL; rootCount: NAT; MM: REAL _ cone.M*cone.M; MMh: REAL _ MM*cone.h; p: Point3d _ localRay.basePt; d: Vector _ localRay.direction; hits _ GetFaceHitFromPool[]; a _ d[1]*d[1] + d[3]*d[3] - MM*d[2]*d[2]; b _ 2*(p[1]*d[1] + p[3]*d[3] + MMh*d[2] - MM*p[2]*d[2]); c _ p[1]*p[1] + p[3]*p[3] - MMh*cone.h + 2*MMh*p[2] - MM*p[2]*p[2]; [roots, rootCount] _ PositiveRoots[a, b, c]; hits.count _ 0; FOR i: NAT IN[1..rootCount] DO y _ p[2] + roots[i]*d[2]; IF cone.yLo <= y AND y<= cone.yHi THEN { vectorT, vectorS: Vector; x, z: REAL; hits.count _ hits.count + 1; hits.params[hits.count] _ roots[i]; x _ p[1] + roots[i]*d[1]; z _ p[3] + roots[i]*d[3]; vectorT _ [z, 0, -x]; vectorS _ [-x, cone.h-y, -z]; hits.normals[hits.count] _ SVVector3d.CrossProduct[vectorT, vectorS]; IF cone.normalPointsOut THEN { IF NOT cone.noseIsUp THEN hits.normals[hits.count] _ SVVector3d.Negate[hits.normals[hits.count]]} ELSE { IF cone.noseIsUp THEN hits.normals[hits.count] _ SVVector3d.Negate[hits.normals[hits.count]]}; }; -- end of root hits ENDLOOP; }; -- end of ConeCast DiskRingCast: PUBLIC PROC [diskRing: DiskRing, localRay: Ray] RETURNS [hits: FaceHitRec] = { x, z, r2, t: REAL; p: Point3d _ localRay.basePt; d: Vector _ localRay.direction; hits _ GetFaceHitFromPool[]; hits.count _ 0; IF Abs[d[2]] < almostZero THEN RETURN; -- ray parallel to plane of disk t _ (diskRing.yPlane - p[2])/d[2]; IF t > 0 THEN { x _ p[1] + t*d[1]; z _ p[3] + t*d[3]; r2 _ x*x + z*z; IF diskRing.rInSquared <= r2 AND r2 <= diskRing.rOutSquared THEN { hits.params[1] _ t; hits.count _ 1; hits.normals[1] _ diskRing.normal; }; }; }; -- end of DiskRingCast CylinderCast: PUBLIC PROC [cylinder: Cylinder, localRay: Ray] RETURNS [hits: FaceHitRec] = { x, y, z: REAL; a, b, c: REAL; p: Point3d _ localRay.basePt; d: Vector _ localRay.direction; roots: ARRAY[1..2] OF REAL; rootCount: NAT; hits _ GetFaceHitFromPool[]; hits.count _ 0; a _ d[1]*d[1] + d[3]*d[3]; b _ 2*(p[1]*d[1] + p[3]*d[3]); c _ p[1]*p[1] + p[3]*p[3] - cylinder.rSquared; [roots, rootCount] _ PositiveRoots[a, b, c]; FOR i: NAT IN[1..rootCount] DO y _ p[2] + roots[i]*d[2]; IF cylinder.yLo <= y AND y <= cylinder.yHi THEN { hits.count _ hits.count + 1; hits.params[hits.count] _ roots[i]; x _ p[1] + roots[i]*d[1]; z _ p[3] + roots[i]*d[3]; hits.normals[hits.count] _ [x/cylinder.r, 0, z/cylinder.r]; IF NOT cylinder.normalPointsOut THEN hits.normals[hits.count] _ SVVector3d.Negate[hits.normals[hits.count]]; }; -- end of y in bounds ENDLOOP; }; -- end of CylinderCast EdgeOnRectCast: PUBLIC PROC [rect: EdgeOnRect, localRay: Ray] RETURNS [hits: FaceHitRec] = { x, y, z, t, denom: REAL; hits _ GetFaceHitFromPool[]; denom _ rect.A*localRay.direction[1] + rect.B*localRay.direction[2]; IF Abs[denom] < almostZero THEN { hits.count _ 0; RETURN}; t _ -(rect.A*localRay.basePt[1] + rect.B*localRay.basePt[2] + rect.D)/denom; IF t <= 0 THEN {hits.count _ 0; RETURN}; x _ localRay.basePt[1] + t*localRay.direction[1]; y _ localRay.basePt[2] + t*localRay.direction[2]; z _ localRay.basePt[3] + t*localRay.direction[3]; IF z > rect.frontZ OR z < rect.backZ THEN {hits.count _ 0; RETURN}; IF rect.valIsY THEN { IF y > rect.valHi OR y < rect.valLo THEN {hits.count _ 0; RETURN}} ELSE { IF x > rect.valHi OR x < rect.valLo THEN {hits.count _ 0; RETURN}}; hits.count _ 1; hits.params[1] _ t; hits.normals[1] _ [rect.A, rect.B, 0]; }; -- end of EdgeOnRectCast globalFaceHitCount: NAT; globalFaceHitMaxCount: NAT = 20; globalFaceHits: ARRAY [1..globalFaceHitMaxCount] OF FaceHitRec; GetFaceHitFromPool: PUBLIC PROC RETURNS [faceHit: FaceHitRec] = { faceHit _ globalFaceHits[globalFaceHitCount]; globalFaceHitCount _ globalFaceHitCount - 1; }; ReturnFaceHitToPool: PUBLIC PROC [faceHit: FaceHitRec] = { globalFaceHitCount _ globalFaceHitCount + 1; globalFaceHits[globalFaceHitCount] _ faceHit; }; Init: PROC = { globalFaceHitCount _ globalFaceHitMaxCount; FOR i: NAT IN[1..globalFaceHitMaxCount] DO globalFaceHits[i] _ NEW[FaceHitRecObj]; ENDLOOP; }; Init[]; END. ªFile: SVFacesCastImpl.mesa Author: Eric Bier on July 4, 1983 10:51 pm Last edited by Bier on July 4, 1983 11:02 pm Contents: Routines for finding the intersections of rays (see CSG.mesa) with simple faces (see SVFaces.mesa) Use only the positive roots. Do one ray-cone intersection and then test for limits. Find the normal for each hit. The cone equations are: x^2 + z^2 = r^2. r = M(h-y). which expands to: x^2 + z^2 = M^2*h^2 - 2*M*h*y + M^2*y^2 Substituting in the ray equations: x = p1 + t*d1; y = p2 + t*d2; z = p3 + t*d3; we have: t^2*(d1d1 + d3d3 - d2d2) + t*(2*p1d1 + 2*p3d3 + 2*Mhd2 - 2*p2d2) + (p1p1 + p3p3 - MhMh + 2*Mhp2 - p2p2) = 0. Is this intersection within the bounds of our cone? calculate the normal A vector clockwise around the circle x^2 + z^2 = x0^2 + z0^2 A vector from the intersection to the apex [0, cone.h, 0] If the cone is nose up, this cross product will point out, else in. Do 1 ray-plane intersection and test to see if it is between the two circles Plane is y = diskRing.yPlane; For ray, y = p2 + t*d2; Solving, diskRing.yPlane = p2 + t*d2; (diskRing.yPlane - p2)/d2 = t; Do one ray-cylinder test and check for any intersections which are between yHi and yLo. The cylinder equations are: x^2 + z^2 = r^2 for yLo <= y <= yHi. The ray equations are, as usual: x = p1 + t*d1; y = p2 + t*d2; z = p3 + t*d3; Plugging ray into cylinder, we get: t^2*(d1d1 + d3d3) + t*(2*p1d1 + 2*p3d3) + p1p1 + p3p3 - r*r = 0; Solve for t; See if y is in bounds Find the normal Do one ray-edge on rect test. The plane of an edge on rect has equation: Ax + By + D = 0; Plugging in the ray equations, we have: A(p1 + td1) + B(p2 + td2) + D = 0; t(Ad1 + Bd2) + Ap1 + Bp2 + D = 0. Solving for t: t = -(Ap1 + Bp2 + D)/(Ad1 + Bd2) Now test to see if intersection point is in the bounds of the edge on rectangle. Its z must be between rect.frontZ and rect.backZ and either its x or its y must be between valLo and valHi depending on valIsY. ÊF– "cedar" style˜Iheadšœ™Iprocšœ*™*Lšœ,™,Lšœl™lL˜šÏk ˜ Lšœ ˜ Lšœ˜Lšœ ˜ Lšœ˜Lšœ˜Lšœ˜—L˜šœ˜Lšœ˜Lšœ˜—Lš˜˜Lšœœœ ˜Lšœ œ˜ Lšœ œœ ˜!Lšœ œ˜(Lšœ œœ ˜!Lšœ œ˜(Lšœœœ˜%Lšœœ˜,Lšœ œ˜!Lšœœœ˜Lšœœœ˜Lšœœ˜"L˜Lšœœœ˜%Lšœœ˜0L˜Lšœœ˜L˜š Ïnœœœœœ˜&Lšœœœœ˜"Lšœ˜—L˜šž œ œœœ œœœ œ˜eJ™Lšœœ˜Lšœ=˜=Lšœœœ˜šœ˜šœœ˜šœœœ˜!Lšœ˜—Lšœ˜Lšœ˜L˜—Lšœ ˜—Lšœ˜Jšœ˜—š œžœœœœ˜QLšœU™ULšœ5™5Lšœ9™9LšœO™OLšœI™ILšœ+™+Lšœ œ˜Lšœœ˜Lšœœœœ˜Lšœœ˜Lšœœœœ˜Lšœœœ˜Lšœ˜Lšœ˜Lšœ˜Lšœœ˜)Lšœ*œ˜8Lšœ6œ˜CLšœ,˜,Lšœ˜šœœœ˜Lšœ3™3Lšœ˜šœœœ˜(Lšœ˜Lšœœ˜Lšœ˜Lšœ#˜#Lšœ™Lšœ˜Lšœ˜Lšœ˜Lšœ<™