DIRECTORY Algebra3d, SVRay, SV3d, SVFaces, SVFacesCast, SVRayTypes, SVVector3d; SVFacesCastImpl: CEDAR PROGRAM IMPORTS Algebra3d, SVRay, SVVector3d EXPORTS SVFacesCast = BEGIN Cone: TYPE = SVFaces.Cone; DiskRing: TYPE = SVFaces.DiskRing; Cylinder: TYPE = SVFaces.Cylinder; EdgeOnRect: TYPE = SVFaces.EdgeOnRect; Point3d: TYPE = SV3d.Point3d; Ray: TYPE = SVRayTypes.Ray; Vector3d: TYPE = SV3d.Vector3d; FaceHitRec: TYPE = REF FaceHitRecObj; FaceHitRecObj: TYPE = SVFacesCast.FaceHitRecObj; almostZero: REAL = 1.0e-12; 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; d: Vector3d; [p, d] _ SVRay.GetLocalRay[localRay]; 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 { 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]; hits.normals[hits.count] _ [x, MM*(cone.h-y), z]; IF NOT cone.normalPointsOut 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; d: Vector3d; [p, d] _ SVRay.GetLocalRay[localRay]; 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; d: Vector3d; roots: ARRAY[1..2] OF REAL; rootCount: NAT; [p, d] _ SVRay.GetLocalRay[localRay]; 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, 0, z]; 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; p: Point3d; d: Vector3d; [p, d] _ SVRay.GetLocalRay[localRay]; hits _ GetFaceHitFromPool[]; denom _ rect.A*d[1] + rect.B*d[2]; IF ABS[denom] < almostZero THEN { hits.count _ 0; RETURN}; t _ -(rect.A*p[1] + rect.B*p[2] + rect.D)/denom; IF t <= 0 THEN {hits.count _ 0; RETURN}; x _ p[1] + t*d[1]; y _ p[2] + t*d[2]; z _ p[3] + t*d[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 Copyright c 1984 by Xerox Corporation. All rights reserved. Last edited by Bier on February 18, 1987 6:27:59 pm PST Contents: Routines for finding the intersections of rays (see SVRay.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 This vector will point out. 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. Ê’– "cedar" style˜Iheadšœ™Jšœ Ïmœ1™