File: SVFacesCastImpl.mesa
Copyright © 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)
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] = {
Use only the positive roots.
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] = {
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.
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
Is this intersection within the bounds of our cone?
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];
Calculate the normal
x ← p[1] + roots[i]*d[1];
z ← p[3] + roots[i]*d[3];
This vector will point out.
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] = {
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;
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;
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;
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
See if y is in bounds
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];
Find the normal
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] = {
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)
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];
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.
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;
};
END.