File: SVObjectCastImplA.mesa
Created July 20, 1982 3:36 pm
Copyright © 1984 by Xerox Corporation. All rights reserved.
Last edited by Bier on June 2, 1987 1:10:51 pm PDT
Contents: Ray-Object intersection tests for block, sphere, and cylinder
DIRECTORY
Algebra3d, Feedback, Rope, SV2d, SV3d, SVCastRays, SVMasterObject, SVMasterObjectTypes, SVModelTypes, SVObjectCast, SVRay, SVSceneTypes, SVVector3d;
SVObjectCastImplA: CEDAR PROGRAM
IMPORTS Algebra3d, SVCastRays, SVRay, Feedback, SVVector3d
EXPORTS SVObjectCast =
BEGIN
BlockHitData: TYPE = REF BlockHitDataObj;
BlockHitDataObj: TYPE = SVMasterObject.BlockHitDataObj;
Classification: TYPE = SVSceneTypes.Classification;
Composite: TYPE = SVSceneTypes.Composite;
CSGTree: TYPE = SVSceneTypes.CSGTree;
LightSourceList: TYPE = SVModelTypes.LightSourceList;
Matrix4by4: TYPE = SV3d.Matrix4by4;
Point2d: TYPE = SV2d.Point2d;
Point3d: TYPE = SV3d.Point3d;
PointSetOp: TYPE = SVSceneTypes.PointSetOp;
Primitive: TYPE = SVSceneTypes.Primitive;
Ray: TYPE = SVSceneTypes.Ray;
Surface: TYPE = REF ANY;
Vector3d: TYPE = SV3d.Vector3d;
RectSurface: TYPE = SVObjectCast.RectSurface;
RectType: TYPE = SVObjectCast.RectType;-- {up, down, front, back, left, right};
TubeSurface: TYPE = SVObjectCast.TubeSurface;
DiskSurface: TYPE = SVObjectCast.DiskSurface;
DiskType: TYPE = SVObjectCast.DiskType;-- {top, bottom};
ShellSurface: TYPE = SVObjectCast.ShellSurface;
HitRec: TYPE = REF HitRecObj;
HitRecObj: TYPE = SVObjectCast.HitRecObj;
HitArray: TYPE = REF HitArrayObj;
HitArrayObj: TYPE = SVObjectCast.HitArrayObj;
SurfaceArray: TYPE = SVSceneTypes.SurfaceArray;
ParameterArray: TYPE = SVSceneTypes.ParameterArray;
InOutArray: TYPE = SVSceneTypes.InOutArray;
NormalArray: TYPE = SVSceneTypes.NormalArray;
BlockData: TYPE = SVMasterObjectTypes.BlockData;
globalHitArray: HitArray;
globalHitArray2: HitArray;
Between: PROC [test, bound1, bound2: REAL] RETURNS [BOOL] = {
min, max: REAL;
min ← MIN[bound1, bound2];
max ← MAX[bound1, bound2];
RETURN[min<= test AND test <= max];
};
BlockCast: PUBLIC PROC [localRay: Ray, prim: Primitive, blockRec: BlockData, surfs: SurfaceArray] RETURNS [class: Classification] = {
Do six ray-plane intersection tests.
Front face is at z=hiZ with loX < x < hiX and loY < y < hiY
Ray parameterization is Z = p[3] + t*d[3]; and so on.
We wish to find the value of t for which z = hiZ
Solving: hiZ = p[3] + t*d[3], or t = (hiZ - p[3])/d[3].
There are 3 possibilities:
1) t comes out negative. The front face is on the wrong half ray. Count this as a miss.
2) Direction[3] is 0. The ray is parallel to the face. Record no intersections.
3) t is valid. Calculate x and y. Check for loX < x < hiX and loY < y < hiY. Record a hit if these are true.
Because x(t) = p[1] + t*d[1]. y(t) = p[2] + t*d[2], the inequalities reduce to:
loX < p[1] + ((hiZ - p[3])/d[3])d[1] < hiX
or loX d[3] < p[1]d[3] + hiZ d[1] - p[3]d[1] < hiX d[3] x test for the front face
loY d[3] < p[2]d[3] + hiZ d[2] - p[3]d[2] < hiY d[3] y test for the front face
loX d[3] < p[1]d[3] + loZ d[1] - p[3]d[1] < hiX d[3] x test for the back face
loY d[3] < p[2]d[3] + loZ d[2] - p[3]d[2] < hiY d[3] y test for the back face

loX d[2] < p[1]d[2] + hiY d[1] - p[2]d[1] < hiX d[2] x test for top face
loZ d[2] < p[3]d[2] + hiY d[3] - p[2]d[3] < hiZ d[2] z test for top face
loX d[2] < p[1]d[2] + loY d[1] - p[2]d[1] < hiX d[2] x test for bottom face
loZ d[2] < p[3]d[2] + loY d[3] - p[2]d[3] < hiZ d[2] z test for bottom face

loY d[1] < p[2]d[1] + hiX d[2] - p[1]d[2] < hiY d[1] y test for right face
loZ d[1] < p[3]d[1] + hiX d[3] - p[1]d[3] < hiZ d[1] z test for right face
loY d[1] < p[2]d[1] + loX d[2] - p[1]d[2] < hiY d[1] y test for left face
loZ d[1] < p[3]d[1] + loX d[3] - p[1]d[3] < hiZ d[1] z test for left face
change the sign of the limit d[*] in d[*] < * < d[*] if d[*] is negative
notice that with the six p(*)d(*) cross terms and the 18 [loX, hiX, loY, hiY, loZ, hiZ]d(*) cross terms, we can do all of these tests
t: REAL;
xd, yd, zd, D: REAL;
hits: HitArray ← globalHitArray;
loX, loY, loZ, hiX, hiY, hiZ: REAL;
p1d2, p1d3, p2d1, p2d3, p3d1, p3d2: REAL;
loZd1, hiZd1, loZd2, hiZd2, loZd3, hiZd3: REAL;
loYd1, hiYd1, loYd2, hiYd2, loYd3, hiYd3: REAL;
loXd1, hiXd1, loXd2, hiXd2, loXd3, hiXd3: REAL;
p: Point3d;
d: Vector3d;
hitCount: NAT ← 0;
doesHit: BOOL;
upS, downS, frontS, backS, rightS, leftS: Surface;
almostZero: REAL ← 1.0e-12;
loX ← blockRec.box.loX; loY ← blockRec.box.loY; loZ ← blockRec.box.loZ;
hiX ← blockRec.box.hiX; hiY ← blockRec.box.hiY; hiZ ← blockRec.box.hiZ;
upS ← surfs[1]; downS ← surfs[2];
leftS ← surfs[3]; backS ← surfs[4]; rightS ← surfs[5]; frontS ← surfs[6];
class ← SVCastRays.GetClassFromPool[];
[p, d] ← SVRay.GetLocalRay[localRay];
p1d2 ← p[1]*d[2]; p1d3 ← p[1]*d[3]; p2d1 ← p[2]*d[1];
p2d3 ← p[2]*d[3]; p3d1 ← p[3]*d[1]; p3d2 ← p[3]*d[2];
loZd1 ← loZ*d[1]; hiZd1 ← hiZ*d[1]; loZd2 ← loZ*d[2]; hiZd2 ← hiZ*d[2]; loZd3 ← loZ*d[3]; hiZd3 ← hiZ*d[3];
loYd1 ← loY*d[1]; hiYd1 ← hiY*d[1]; loYd2 ← loY*d[2]; hiYd2 ← hiY*d[2]; loYd3 ← loY*d[3]; hiYd3 ← hiY*d[3];
loXd1 ← loX*d[1]; hiXd1 ← hiX*d[1]; loXd2 ← loX*d[2]; hiXd2 ← hiX*d[2]; loXd3 ← loX*d[3]; hiXd3 ← hiX*d[3];
front face and back face
D ← Abs[d[3]];
IF D > almostZero THEN {
front face
xd ← p1d3 + hiZd1 - p3d1;
IF Between[xd, loXd3, hiXd3] THEN {
yd ← p2d3 + hiZd2 - p3d2;
IF Between[yd, loYd3, hiYd3]
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit THEN {
t ← (hiZ - p[3])/d[3];
IF t>0 THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, frontS, [0,0,1]];
};
};
back face
xd ← p1d3 + loZd1 - p3d1;
IF Between[xd, loXd3, hiXd3] THEN {
yd ← p2d3 + loZd2 - p3d2;
IF Between[yd, loYd3, hiYd3]
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit THEN {
t ← (loZ - p[3])/d[3];
IF t>0 THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, backS, [0,0,-1]];
};
};
};
up face and down face
D ← Abs[d[2]];
IF D > almostZero THEN {
up face
xd ← p1d2 + hiYd1 - p2d1;
IF Between[xd, loXd2, hiXd2] THEN {
zd ← p3d2 + hiYd3 - p2d3;
IF Between[zd, loZd2, hiZd2]
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit THEN {
t ← (hiY - p[2])/d[2];
IF t>0 THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, upS, [0,1,0]];
};
};
down face
xd ← p1d2 + loYd1 - p2d1;
IF Between[xd, loXd2, hiXd2] THEN {
zd ← p3d2 + loYd3 - p2d3;
IF Between[zd, loZd2, hiZd2]
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit THEN {
t ← (loY - p[2])/d[2];
IF t>0 THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, downS, [0,-1,0]];
};
};
};
left face and right face
D ← Abs[d[1]];
IF D > almostZero THEN {
right face
yd ← p2d1 + hiXd2 - p1d2;
IF Between[yd, loYd1, hiYd1] THEN {
zd ← p3d1 + hiXd3 - p1d3;
IF Between[zd, loZd1, hiZd1]
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit THEN {
t ← (hiX - p[1])/d[1];
IF t>0 THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, rightS, [1,0,0]];
};
};
left face
yd ← p2d1 + loXd2 - p1d2;
IF Between[yd, loYd1, hiYd1] THEN {
zd ← p3d1 + loXd3 - p1d3;
IF Between[zd, loZd1, hiZd1]
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit THEN {
t ← (loX - p[1])/d[1];
IF t>0 THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, leftS, [-1,0,0]];
};
};
};
SELECT hitCount FROM
0 => { -- ray misses block completely
class.count ← 0;
class.classifs[1] ← FALSE;
};
1 => { -- ray skims an edge and for some reason is noticed by only one surface or ray originates from inside the cube. If the later is true, class has one hit and classifs are TRUE, FALSE.
IF loX < p[1] AND p[1] < hiX AND loY < p[2] AND p[2] < hiY AND loZ < p[3] AND p[3] < hiZ THEN { -- ray starts from inside the cube.
class.count ← 1;
class.params[1] ← hits[1].t;
class.surfaces[1] ← hits[1].surf;
class.stills[1] ← NARROW[hits[1].surf, BlockHitData].still;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.normals[1] ← hits[1].normal;
class.primitives[1] ← prim;
}
ELSE { -- Count as a miss.
SVCastRays.MakeClassAMiss[class];
};
};
2 => { -- ray enters and exits block properly
The hit with the lower value of t is the entry point. The other is the exit.
class.count ← 2;
IF hits[1].t < hits[2].t THEN {
class.params[1] ← hits[1].t;
class.params[2] ← hits[2].t;
class.surfaces[1] ← hits[1].surf;
class.surfaces[2] ← hits[2].surf;
class.stills[1] ← NARROW[hits[1].surf, BlockHitData].still;
class.stills[2] ← NARROW[hits[2].surf, BlockHitData].still;
class.normals[1] ← hits[1].normal;
class.normals[2] ← hits[2].normal;}
ELSE {
class.params[2] ← hits[1].t;
class.params[1] ← hits[2].t;
class.surfaces[2] ← hits[1].surf;
class.surfaces[1] ← hits[2].surf;
class.stills[2] ← NARROW[hits[1].surf, BlockHitData].still;
class.stills[1] ← NARROW[hits[2].surf, BlockHitData].still;
class.normals[2] ← hits[1].normal;
class.normals[1] ← hits[2].normal;};
class.classifs[1] ← FALSE; class.classifs[2] ← TRUE; class.classifs[3] ← FALSE;
class.primitives[1] ← class.primitives[2] ← prim;
};
3 => { -- ray hits a vertex and is noticed by all 3 surfaces. Treat as two hits for now using the frontmost and rearmost values of class.params.
perm: Permutation ← Rank[[hits[1].t, hits[2].t, hits[3].t]];
class.count ← 2;
class.params[1] ← hits[perm[1]].t;
class.params[2] ← hits[perm[3]].t;
class.surfaces[1] ← hits[perm[1]].surf;
class.surfaces[2] ← hits[perm[3]].surf;
class.stills[1] ← NARROW[hits[1].surf, BlockHitData].still;
class.stills[2] ← NARROW[hits[3].surf, BlockHitData].still;
class.normals[1] ← hits[perm[1]].normal;
class.normals[2] ← hits[perm[3]].normal;
class.classifs[1] ← FALSE; class.classifs[2] ← TRUE; class.classifs[3] ← FALSE;
class.primitives[1] ← class.primitives[2] ← prim;
Feedback.AppendTypescriptRaw[$Solidviews, "Point of interest: Some block was hit 3 times", oneLiner];
};
ENDCASE => SIGNAL WierdHitCount;
}; -- end of BlockCast
Here is the version that works for cubes of side 2 centered on the origin.
BlockCast: PUBLIC PROC [localRay: Ray, prim: Primitive, blockRec: BlockData, surfs: SurfaceArray] RETURNS [class: Classification] = {
Do six ray-plane intersection tests.
Front face is at z=1.0 with -1.0 < x < 1.0 and -1.0 < y < 1.0
Ray parameterization is Z = p[3] + t*d[3]; and so on.
We wish to find the value of t for which z = 1.0
Solving: 1.0 = p[3] + t*d[3], or t = (1 - p[3])/d[3].
There are 3 possibilities:
1) t comes out negative. The front face is on the wrong half ray. Count this as a miss.
2) Direction[3] is 0. The ray is parallel to the face. Record no intersections.
3) t is valid. Calculate x and y. Check for -1.0 < x < 1.0 and -1.0 < y < 1.0. Record a hit if these are true.
Because x(t) = p[1] + t*d[1]. y(t) = p[2] + t*d[2], the inequalities reduce to:
-1 < p[1] + ((1 - p[3])/d[3])d[1] < 1
or -d[3] < p[1]d[3] + d[1] - p[3]d[1] < d[3] x test for the front face
-d[3] < p[2]d[3] + d[2] - p[3]d[2] < d[3] y test for the front face
-d[3] < p[1]d[3] - d[1] - p[3]d[1] < d[3] x test for the back face
-d[3] < p[2]d[3] - d[2] - p[3]d[2] < d[3] y test for the back face

-d[2] < p[1]d[2] + d[1] - p[2]d[1] < d[2] x test for top face
-d[2] < p[3]d[2] + d[3] - p[2]d[3] < d[2] z test for top face
-d[2] < p[1]d[2] - d[1] - p[2]d[1] < d[2] x test for bottom face
-d[2] < p[3]d[2] - d[3] - p[2]d[3] < d[2] z test for bottom face

-d[1] < p[2]d[1] + d[2] - p[1]d[2] < d[1] y test for right face
-d[1] < p[3]d[1] + d[3] - p[1]d[3] < d[1] z test for right face
-d[1] < p[2]d[1] - d[2] - p[1]d[2] < d[1] y test for right face
-d[1] < p[3]d[1] - d[3] - p[1]d[3] < d[1] z test for right face
change the sign of the limit d[*] in d[*] < * < d[*] if d[*] is negative
notice that with the six p(*)d(*) cross terms, we can do all of these tests
t: REAL;
xd, yd, zd, D: REAL;
hits: HitArray ← globalHitArray;
p1d2, p1d3, p2d1, p2d3, p3d1, p3d2: REAL;
p: Point3d;
d: Vector3d;
hitCount: NAT ← 0;
doesHit: BOOL;
upS, downS, frontS, backS, rightS, leftS: Surface;
almostZero: REAL ← 1.0e-12;
upS ← surfs[1]; downS ← surfs[2];
leftS ← surfs[3]; backS ← surfs[4]; rightS ← surfs[5]; frontS ← surfs[6];
class ← SVCastRays.GetClassFromPool[];
[p, d] ← SVRay.GetLocalRay[localRay];
p1d2 ← p[1]*d[2];
p1d3 ← p[1]*d[3];
p2d1 ← p[2]*d[1];
p2d3 ← p[2]*d[3];
p3d1 ← p[3]*d[1];
p3d2 ← p[3]*d[2];
front face and back face
D ← Abs[d[3]];
IF D > almostZero THEN {
front face
xd ← p1d3 + d[1] - p3d1;
IF -D <= xd AND xd <= D THEN {
yd ← p2d3 + d[2] - p3d2;
IF -D <= yd AND yd <= D
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit THEN {
t ← (1.0 - p[3])/d[3];
IF t>0 THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, frontS, [0,0,1]];
};
};
back face
xd ← p1d3 - d[1] - p3d1;
IF -D <= xd AND xd <= D THEN {
yd ← p2d3 - d[2] - p3d2;
IF -D <= yd AND yd <= D
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit THEN {
t ← (-1.0 - p[3])/d[3];
IF t>0 THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, backS, [0,0,-1]];
};
};
};
up face and down face
D ← Abs[d[2]];
IF D > almostZero THEN {
up face
xd ← p1d2 + d[1] - p2d1;
IF -D <= xd AND xd <= D THEN {
zd ← p3d2 + d[3] - p2d3;
IF -D <= zd AND zd <= D
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit THEN {
t ← (1.0 - p[2])/d[2];
IF t>0 THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, upS, [0,1,0]];
};
};
down face
xd ← p1d2 - d[1] - p2d1;
IF -D <= xd AND xd <= D THEN {
zd ← p3d2 - d[3] - p2d3;
IF -D <= zd AND zd <= D
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit THEN {
t ← (-1.0 - p[2])/d[2];
IF t>0 THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, downS, [0,-1,0]];
};
};
};
left face and right face
D ← Abs[d[1]];
IF D > almostZero THEN {
right face
yd ← p2d1 + d[2] - p1d2;
IF -D <= yd AND yd <= D THEN {
zd ← p3d1 + d[3] - p1d3;
IF -D <= zd AND zd <= D
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit THEN {
t ← (1.0 - p[1])/d[1];
IF t>0 THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, rightS, [1,0,0]];
};
};
left face
yd ← p2d1 - d[2] - p1d2;
IF -D <= yd AND yd <= D THEN {
zd ← p3d1 - d[3] - p1d3;
IF -D <= zd AND zd <= D
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit THEN {
t ← (-1.0 - p[1])/d[1];
IF t>0 THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, leftS, [-1,0,0]];
};
};
};
SELECT hitCount FROM
0 => { -- ray misses block completely
class.count ← 0;
class.classifs[1] ← FALSE;
};
1 => { -- ray skims an edge and for some reason is noticed by only one surface or ray originates from inside the cube. If the later is true, class has one hit and classifs are TRUE, FALSE.
IF p[1] > -1.0 AND p[1] < 1.0 AND p[2] > -1.0 AND p[2] < 1.0 AND p[3] > -1.0 AND p[3] < 1.0 THEN { -- ray starts from inside the cube.
class.count ← 1;
class.params[1] ← hits[1].t;
class.surfaces[1] ← hits[1].surf;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.normals[1] ← hits[1].normal;
class.primitives[1] ← prim;
}
ELSE { -- Count as a miss.
SVCastRays.MakeClassAMiss[class];
};
};
2 => { -- ray enters and exits block properly
The hit with the lower value of t is the entry point. The other is the exit.
class.count ← 2;
IF hits[1].t < hits[2].t THEN {
class.params[1] ← hits[1].t;
class.params[2] ← hits[2].t;
class.surfaces[1] ← hits[1].surf;
class.surfaces[2] ← hits[2].surf;
class.normals[1] ← hits[1].normal;
class.normals[2] ← hits[2].normal;}
ELSE {
class.params[2] ← hits[1].t;
class.params[1] ← hits[2].t;
class.surfaces[2] ← hits[1].surf;
class.surfaces[1] ← hits[2].surf;
class.normals[2] ← hits[1].normal;
class.normals[1] ← hits[2].normal;};
class.classifs[1] ← FALSE; class.classifs[2] ← TRUE; class.classifs[3] ← FALSE;
class.primitives[1] ← class.primitives[2] ← prim;
};
3 => { -- ray hits a vertex and is noticed by all 3 surfaces. Treat as two hits for now using the frontmost and rearmost values of class.params.
perm: Permutation ← Rank[[hits[1].t, hits[2].t, hits[3].t]];
class.count ← 2;
class.params[1] ← hits[perm[1]].t;
class.params[2] ← hits[perm[3]].t;
class.surfaces[1] ← hits[perm[1]].surf;
class.surfaces[2] ← hits[perm[3]].surf;
class.normals[1] ← hits[perm[1]].normal;
class.normals[2] ← hits[perm[3]].normal;
class.classifs[1] ← FALSE; class.classifs[2] ← TRUE; class.classifs[3] ← FALSE;
class.primitives[1] ← class.primitives[2] ← prim;
Feedback.AppendTypescriptRaw[$Solidviews, "Point of interest: Some block was hit 3 times", oneLiner];
};
ENDCASE => SIGNAL WierdHitCount;
}; -- end of BlockCast
WierdHitCount: SIGNAL = CODE;
Permutation: TYPE = ARRAY [1..3] OF NAT;
Rank: PROC [t: ARRAY[1..3] OF REAL] RETURNS [perm: Permutation] = {
tempT: REAL;
tempP: NAT;
perm ← [1, 2, 3];
FOR i: NAT IN[1..2] DO
FOR j: NAT IN[1..3-i] DO
IF t[j] > t[j+1] THEN {
tempT ← t[j];
tempP ← perm[j];
t[j] ← t[j+1];
perm[j] ← perm[j+1];
t[j+1] ← tempT;
perm[j+1] ← tempP;};
ENDLOOP;
ENDLOOP;
}; -- end of Rank
PositiveRoots: PROCEDURE [a, b, c: REAL] RETURNS [rootArray: ARRAY[1..2] OF REAL, rootCount, totalRealRoots: NAT] = {
Use only the positive roots. Preserves the order of the roots, so they will still be in increasing order.
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;
};
AllRoots: PROCEDURE [a, b, c: REAL] RETURNS [rootArray: ARRAY[1..2] OF REAL, positiveRootCount, totalRealRoots: NAT] = {
Count the positive roots. Preserves the order of the roots, so they will still be in increasing order.
positiveRootCount ← 0;
[rootArray, totalRealRoots] ← Algebra3d.QuadraticFormula[a, b, c];
IF totalRealRoots = 0 THEN RETURN;
FOR i: NAT IN [1..totalRealRoots] DO
IF rootArray[i] <= 0 THEN positiveRootCount ← positiveRootCount + 1;
ENDLOOP;
};
SphereCast: PUBLIC PROC [localRay: Ray, thisShell: Surface, prim: Primitive] RETURNS [class: Classification] = {
Do one ray-quadric test.
sphere has radius = 1.0; Equation: x^2 + y^2 + z^2 = r^2.
Find 1 or 2 values of t which satisfy this equation.
r: REAL = 1.0;
The ray has the form:
x = p1 + t*d1
y = p2 + t*d2
z = p3 + t*d3
hence, at the intersection of ray and sphere, t satisfies:
(p1 + t*d1)^2 + (p2 + t*d2)^2 + (p3 + t*d3)^2 = r^2
which reduces to:
t^2 * (d1^2 + d2^2 + d3^2) + t * 2(p1*d1 + p2*d2 + p3*d3) +
(p1^2 + p2^2 + p3^2 - r^2) = 0;
Let p = [p1, p2, p3], and d = [d1, d2, d3] then in dot product notation:
t^2 * d.d + t * 2 * d.p + p.p - r^2 = 0. Use quadratic formula
p: Point3d;
d: Vector3d;
tArray: ARRAY[1..2] OF REAL;
rootCount: NAT;
dp: REAL;
dp2: REAL;
dd: REAL;
pp: REAL;
t1, t2: REAL;
[p, d] ← SVRay.GetLocalRay[localRay];
dp ← SVVector3d.DotProduct[d,p];
dp2 ← dp*dp;
dd ← SVVector3d.DotProduct[d,d];
pp ← SVVector3d.DotProduct[p,p];
class ← SVCastRays.GetClassFromPool[];
[tArray, rootCount, ----] ← PositiveRoots[dd,2*dp,pp-r*r];
SELECT rootCount FROM
0 => {
class.count ← 0;
class.classifs[1] ← FALSE;
};
1 => {
If the ray originates inside the sphere, treat as a single hit. Otherwise treat as a miss.
IF pp <1.0 THEN { -- ray originates inside the sphere
class.count ← 1;
t1 ← tArray[1];
class.classifs[1] ← TRUE;
class.classifs[2] ← FALSE;
class.params[1] ← t1;
class.surfaces[1] ← thisShell;
class.normals[1][1] ← p[1] + d[1]*t1;
class.normals[1][2] ← p[2] + d[2]*t1;
class.normals[1][3] ← p[3] + d[3]*t1;
class.primitives[1] ← prim;
}
ELSE SVCastRays.MakeClassAMiss[class];
};
2 => {
class.count ← 2;
t1 ← tArray[1]; t2 ← tArray[2]; -- since tArray[1] < tArray[2] is guaranteed
class.params[1] ← t1; class.params[2] ← t2;
class.surfaces[1] ← thisShell; class.surfaces[2] ← thisShell;
class.classifs[1] ← FALSE; class.classifs[2] ← TRUE; class.classifs[3] ← FALSE;
class.primitives[1] ← class.primitives[2] ← prim;
class.normals[1][1] ← p[1] + d[1]*t1;
class.normals[1][2] ← p[2] + d[2]*t1;
class.normals[1][3] ← p[3] + d[3]*t1;
class.normals[2][1] ← p[1] + d[1]*t2;
class.normals[2][2] ← p[2] + d[2]*t2;
class.normals[2][3] ← p[3] + d[3]*t2;
};
ENDCASE => ERROR;
};
CylinderCast: PUBLIC PROC [localRay: Ray, surfaces: SurfaceArray, prim: Primitive] RETURNS [class: Classification] = {
2 ray-plane tests and 1 ray-tube test
Cylinder Equation: x^2 + z^2 = r^2 for yLow < y < yHigh
As usual, ray equation is:x = x0 + t*Dx
y = y0 + t*Dy
z = z0 + t*Dz
At intersection, t satisfies:
(x0 + tDx)^2 + (z0 + tDz)^2 = r^2 for yLow <= y <= yHigh
Which reduces to:
t^2 * (Dx^2 + Dz^2) + t * 2 * (x0*Dx + z0 * Dz) + (x0^2 + z0^2 - r^2) = 0;
We find the t's which hit this infinite cylinder and then check for the y bounds:
yLow <= (y0 + t*Dy) <= yHigh
yLow: REAL = -1.0; yHigh: REAL = 1.0;
r: REAL = 1.0;(built into code)
p: Point3d;
d: Vector3d;
tArray: ARRAY[1..2] OF REAL;
a, b, c: REAL;-- the quadratic formula variables
positiveRootCount, totalRealRoots: NAT;
almostZero: REAL ← 1.0e-15;
t, fromYaxis2, y, tempT: REAL;
tubeHits: HitArray ← globalHitArray;
tubeHitCount: NAT ← 0;
planeHits: HitArray ← globalHitArray2;
planeHitCount: NAT ← 0;
totalHits: NAT;
topS, bottomS, tubeS, tempSurface: Surface;
thisShell: Surface;
tubeNormal, tempNormal: Vector3d;
p1d2, p2d1, p2d3, p3d2, d2d2, xd2, zd2: REAL;
[p, d] ← SVRay.GetLocalRay[localRay];
topS ← surfaces[1]; bottomS ← surfaces[2]; tubeS ← surfaces[3];
class ← SVCastRays.GetClassFromPool[];
Do tube test.
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] - 1.0;-- r*r = 1.0
[tArray, positiveRootCount, totalRealRoots] ← AllRoots[a, b, c];
Do 2 ray-plane tests
.top face and bottom face
IF Abs[d[2]] > almostZero THEN {
top face.
ray hits plane of top when t = (yHigh - p[2])/d[2]; yHigh = 1.
at which point, x = p[1] + t*d[1] and z ← p[3] + t*d[3]
so x = p[1] + (1 - p[2])d[1]/d[2] and z = p[3] + (1- p[2])d[3]/d[2] at the top face
and x = p[1] + (-1 - p[2])d[1]/d[2] and z = p[3] + (-1- p[2])d[3]/d[2] at the bottom face
therefore x*d[2] = p1d2 + d1 - d1p2 and z*d2 = p3d2 + d3 - d3p2 at the top face
and x*d[2] = p1d2 - d1 - d1p2 and z*d2 = p3d2 - d3 - d3p2 at the bottom face
the point is in the disk if x^2 + z^2 <= r^2 (=1)
or if (xd2)^2 + (zd2)^2 <= (d2)^2
This methods requires 5 mults instead of the 2 divides and 8 mults for the naive method which calculates t first. Add 1 divide to calculate t if point is found to be in disk.
p1d2 ← p[1]*d[2];
p2d1 ← p[2]*d[1];
p2d3 ← p[2]*d[3];
p3d2 ← p[3]*d[2];
d2d2 ← d[2]*d[2];
xd2 ← p1d2 + d[1] - p2d1;
zd2 ← p3d2 + d[3] - p2d3;
fromYaxis2 ← xd2*xd2 + zd2*zd2;
IF fromYaxis2 < d2d2-- r*r = 1.0
THEN {
t ← (1 - p[2])/d[2];
IF t > 0 THEN {
planeHitCount ← planeHitCount + 1;
planeHits[planeHitCount]^ ← [t, topS, [0,1.0,0] ];
};
};
bottom face.
xd2 ← p1d2 - d[1] - p2d1;
zd2 ← p3d2 - d[3] - p2d3;
fromYaxis2 ← xd2*xd2 + zd2*zd2;
IF fromYaxis2 < d2d2-- r*r = 1.0
THEN {
t ← (yLow - p[2])/d[2];
IF t > 0 THEN {
planeHitCount ← planeHitCount + 1;
planeHits[planeHitCount]^ ← [t, bottomS, [0,-1.0,0] ]
};
}; -- end check ray-plane intersections
};
Sort the plane hits.
IF planeHitCount = 2 THEN {
IF planeHits[1].t > planeHits[2].t THEN {
tempT ← planeHits[1].t;
tempSurface ← planeHits[1].surf;
tempNormal ← planeHits[1].normal;
planeHits[1].t ← planeHits[2].t;
planeHits[1].surf ← planeHits[2].surf;
planeHits[1].normal ← planeHits[2].normal;
planeHits[2].t ← tempT;
planeHits[2].surf ← tempSurface;
planeHits[2].normal ← tempNormal;
};
};
totalHits ← planeHitCount + positiveRootCount;
SELECT totalRealRoots FROM
0 => {-- ray does not hit tube but may pass through the planes on the ends
SELECT planeHitCount FROM
0=> {
class.count ← 0;
class.classifs[1] ← FALSE;
};
1=> {
Perhaps a ray is tangent to the tube at the top or bottom of the cylinder. More likely, the ray originates inside of the cylinder. See if this is the case.
IF p[2] > yLow AND p[2] < yHigh AND c < 0.0 THEN {
class.count ← 1;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.params[1] ← planeHits[1].t;
class.normals[1] ← planeHits[1].normal;
class.surfaces[1] ← planeHits[1].surf;
class.primitives[1] ← prim;
}
ELSE SVCastRays.MakeClassAMiss[class];
};
2 => {
class.count ← 2;
class.params[1] ← planeHits[1].t; -- we know planeHits[1].t < planeHits[2].t
class.params[2] ← planeHits[2].t;
class.surfaces[1] ← planeHits[1].surf;
class.surfaces[2] ← planeHits[2].surf;
class.normals[1] ← planeHits[1].normal;
class.normals[2] ← planeHits[2].normal;
class.classifs[1] ← FALSE; class.classifs[2] ← TRUE; class.classifs[3] ← FALSE;
class.primitives[1] ← class.primitives[2] ← prim;
}; -- end 0 tube, 2 disk
ENDCASE => ERROR;
}; -- end of zero tube hits
1 => {
Ray just scratches tube. Count as a miss.
SVCastRays.MakeClassAMiss[class];
y ← p[2] + tArray[1]*d[2];
IF yLow <= y AND y<= yHigh THEN { -- in bounds surface scratch
IF planeHitCount > 0 THEN ERROR; -- could actually happen if d[1]=d[3] = 0;
class.count ← 1;
class.classifs[1] ← FALSE;
class.classifs[2] ← FALSE;
class.params[1] ← tArray[1];
class.surfaces[1] ← thisShell;
class.primitives[1] ← prim;
Normal has x and z components equal to point on the cylinder. y part is 0;
class.normals[1][1] ← p[1] + tArray[1]*d[1];
class.normals[1][2] ← 0;
class.normals[1][3] ← p[3] + tArray[1]*d[3];
}
ELSE{ -- no hit
class.classifs[1] ← FALSE;
};
};
2 => {
Several cases. 0,1 or 2 of these roots may be beyond the end of the cylinder (or on the wrong side of the ray).
If both inside, plane intersections are an error.
If one inside, planeCount should be one.
If none inside, planeCount should be two (if ray hits solid) or zero (if it doesn't).
Calculate number of tube intersections inside cylinder.
tubeHitCount ← 0;
y ← p[2] + tArray[1]*d[2];
IF yLow <= y AND y <= yHigh THEN {-- first tube hit in bounds
tubeHitCount ← tubeHitCount + 1;
tubeNormal[1] ← p[1] + tArray[1]*d[1];
tubeNormal[2] ← 0;
tubeNormal[3] ← p[3] + tArray[1]*d[3];
tubeHits[tubeHitCount]^ ← [tArray[1], tubeS, tubeNormal]};
y ← p[2] + tArray[2]*d[2];
IF yLow <= y AND y <= yHigh THEN {-- second tube hit in bounds
tubeHitCount ← tubeHitCount + 1;
tubeNormal[1] ← p[1] + tArray[2]*d[1];
tubeNormal[2] ← 0;
tubeNormal[3] ← p[3] + tArray[2]*d[3];
tubeHits[tubeHitCount]^ ← [tArray[2], tubeS, tubeNormal]};
SELECT tubeHitCount FROM
2 => { -- both intersections are at the tube
IF planeHitCount >1 THEN ERROR; -- planeHitCount could be 1 if the ray enters in the middle of one tube and exits at a brim. Ignore the plane intersection.
At this point, if any of the intersections are on the wrong side of the ray, assume the ray starts inside:
IF tArray[1] <= 0 THEN {
IF tArray[2] <= 0 THEN SVCastRays.MakeClassAMiss[class]
ELSE {
class.count ← 1;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.params[1] ← tubeHits[2].t;
class.normals[1] ← tubeHits[2].normal;
class.surfaces[1] ← tubeS;
class.primitives[1] ← prim;
};
RETURN;
};
class.count ← 2;
class.params[1] ← tubeHits[1].t; -- we know that tubeHits[1].t < tubeHits[2].t
class.params[2] ← tubeHits[2].t;
class.normals[1] ← tubeHits[1].normal;
class.normals[2] ← tubeHits[2].normal;
class.surfaces[1] ← tubeS; class.surfaces[2] ← tubeS;
class.classifs[1] ← FALSE; class.classifs[2] ← TRUE; class.classifs[3] ← FALSE;
class.primitives[1] ← class.primitives[2] ← prim;
};
1 => { -- one line intersection on the tube. The other is out of bounds. This must be an angled ray. It should either originate from inside or hit a disk.
SELECT planeHitCount FROM
0 => { -- originate from inside?
IF p[2] > yLow AND p[2] < yHigh AND c < 0.0 THEN {
class.count ← 1;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.params[1] ← planeHits[1].t;
class.normals[1] ← planeHits[1].normal;
class.surfaces[1] ← planeHits[1].surf;
class.primitives[1] ← prim;
}
ELSE SVCastRays.MakeClassAMiss[class];
};
1 => {
If either the tubeHit or the planeHit is negative, the ray originates from inside of the cylinder.
IF tubeHits[1].t <= 0 THEN {
IF p[2] > yLow AND p[2] < yHigh AND c < 0.0 THEN {
class.count ← 1;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.params[1] ← planeHits[1].t;
class.normals[1] ← planeHits[1].normal;
class.surfaces[1] ← planeHits[1].surf;
class.primitives[1] ← prim;
}
ELSE SVCastRays.MakeClassAMiss[class];
RETURN;
};
IF planeHits[1].t <=0 THEN {
IF p[2] > yLow AND p[2] < yHigh AND c < 0.0 THEN {
class.count ← 1;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.params[1] ← tubeHits[1].t;
class.normals[1] ← tubeHits[1].normal;
class.surfaces[1] ← tubeS;
class.primitives[1] ← prim;
}
ELSE SVCastRays.MakeClassAMiss[class];
RETURN;
};
class.count ← 2;
IF tubeHits[1].t < planeHits[1].t THEN {
class.params[1] ← tubeHits[1].t;
class.params[2] ← planeHits[1].t;
class.surfaces[1] ← tubeHits[1].surf;
class.surfaces[2] ← planeHits[1].surf;
class.normals[1] ← tubeHits[1].normal;
class.normals[2] ← planeHits[1].normal;}
ELSE {
class.params[2] ← tubeHits[1].t;
class.params[1] ← planeHits[1].t;
class.surfaces[2] ← tubeHits[1].surf;
class.surfaces[1] ← planeHits[1].surf;
class.normals[2] ← tubeHits[1].normal;
class.normals[1] ← planeHits[1].normal;};
class.classifs[1] ← FALSE; class.classifs[2] ← TRUE; class.classifs[3] ← FALSE;
class.primitives[1] ← class.primitives[2] ← prim;
};
2 => {-- we split a disk-plane boundary. Use both plane hits.
If either of the plane hits is on the wrong side of the ray, the ray originates from inside of the cylinder.
IF planeHits[1].t <=0 THEN {
IF p[2] > yLow AND p[2] < yHigh AND c < 0.0 THEN {
class.count ← 1;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.params[1] ← planeHits[2].t;
class.normals[1] ← planeHits[2].normal;
class.surfaces[1] ← planeHits[2].surf;
class.primitives[1] ← prim;
}
ELSE SVCastRays.MakeClassAMiss[class];
RETURN;
};
IF planeHits[2].t <= 0 THEN {
IF p[2] > yLow AND p[2] < yHigh AND c < 0.0 THEN {
class.count ← 1;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.params[1] ← planeHits[1].t;
class.normals[1] ← planeHits[1].normal;
class.surfaces[1] ← planeHits[1].surf;
class.primitives[1] ← prim;
}
ELSE SVCastRays.MakeClassAMiss[class];
RETURN;
};
class.count ← 2;
IF planeHits[1].t < planeHits[2].t THEN {
class.params[1] ← planeHits[1].t;
class.params[2] ← planeHits[2].t;
class.surfaces[1] ← planeHits[1].surf;
class.surfaces[2] ← planeHits[2].surf;
class.normals[1] ← planeHits[1].normal;
class.normals[2] ← planeHits[2].normal;}
ELSE {
class.params[2] ← planeHits[1].t;
class.params[1] ← planeHits[2].t;
class.surfaces[2] ← planeHits[1].surf;
class.surfaces[1] ← planeHits[2].surf;
class.normals[2] ← planeHits[1].normal;
class.normals[1] ← planeHits[2].normal;};
class.classifs[1] ← FALSE; class.classifs[2] ← TRUE; class.classifs[3] ← FALSE;
class.primitives[1] ← class.primitives[2] ← prim;
};
ENDCASE => ERROR;-- can't have four hits total (for now)
};
0 => { -- both intersections are at the disk faces
SELECT planeHitCount FROM
0 => { -- a miss
SVCastRays.MakeClassAMiss[class];
};
1 => { -- perhaps skims a rim but is only noticed by the disk, not the tube. Or perhaps the ray originates from inside the cylinder just below a disk face.
IF p[2] > yLow AND p[2] < yHigh AND c < 0.0 THEN {
class.count ← 1;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.params[1] ← planeHits[1].t;
class.normals[1] ← planeHits[1].normal;
class.surfaces[1] ← planeHits[1].surf;
class.primitives[1] ← prim;
}
ELSE SVCastRays.MakeClassAMiss[class];
};
2 => { -- goes through both disk faces
If either of the plane hits is on the wrong side of the ray, the ray originates from inside of the cylinder. Count as a miss for now.
IF planeHits[1].t <=0 THEN {
IF p[2] > yLow AND p[2] < yHigh AND c < 0.0 THEN {
class.count ← 1;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.params[1] ← planeHits[2].t;
class.normals[1] ← planeHits[2].normal;
class.surfaces[1] ← planeHits[2].surf;
class.primitives[1] ← prim;
}
ELSE SVCastRays.MakeClassAMiss[class];
RETURN;
};
IF planeHits[2].t <= 0 THEN {
IF p[2] > yLow AND p[2] < yHigh AND c < 0.0 THEN {
class.count ← 1;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.params[1] ← planeHits[1].t;
class.normals[1] ← planeHits[1].normal;
class.surfaces[1] ← planeHits[1].surf;
class.primitives[1] ← prim;
}
ELSE SVCastRays.MakeClassAMiss[class];
RETURN;
};
class.count ← 2;
IF planeHits[1].t < planeHits[2].t THEN {
class.params[1] ← planeHits[1].t;
class.params[2] ← planeHits[2].t;
class.surfaces[1] ← planeHits[1].surf;
class.surfaces[2] ← planeHits[2].surf;
class.normals[1] ← planeHits[1].normal;
class.normals[2] ← planeHits[2].normal;}
ELSE {
class.params[2] ← planeHits[1].t;
class.params[1] ← planeHits[2].t;
class.surfaces[2] ← planeHits[1].surf;
class.surfaces[1] ← planeHits[2].surf;
class.normals[2] ← planeHits[1].normal;
class.normals[1] ← planeHits[2].normal;};
class.classifs[1] ← FALSE; class.classifs[2] ← TRUE; class.classifs[3] ← FALSE;
class.primitives[1] ← class.primitives[2] ← prim;
};
ENDCASE => ERROR; -- can't have more than 2 plane hits
End of zero tube hits of 2 tube roots.
};
ENDCASE => ERROR; -- end of select on tubeHitCount with 2 tube roots
}; -- end of 2 tube roots
ENDCASE => ERROR; -- end of select on number of tube roots
}; -- end of procedure CylinderCast
Abs: PROC [r: REAL] RETURNS [REAL] = {
RETURN [IF r>=0 THEN r ELSE -r];
};
Init: PROC = {
Create a global HitArray.
globalHitArray ← NEW[HitArrayObj];
globalHitArray2 ← NEW[HitArrayObj];
FOR i: NAT IN[1..SVObjectCast.globalObjDepth] DO
globalHitArray[i] ← NEW[HitRecObj];
globalHitArray2[i] ← NEW[HitRecObj];
ENDLOOP;
};
Init[];
END.