Point3d: TYPE = Matrix3d.Point3d;
Point2d: TYPE = Matrix3d.Point2d;
Matrix4by4: TYPE = Matrix3d.Matrix4by4;
Vector: TYPE = SVVector3d.Vector;
Surface: TYPE = REF ANY;
LightSourceList: TYPE = DisplayList3d.LightSourceList;
Primitive: TYPE = REF PrimitiveObj;
PrimitiveObj: TYPE = CSG.PrimitiveObj;
RectSurface: TYPE = REF RectSurfaceObj;
RectSurfaceObj: TYPE = ObjectCast.RectSurfaceObj;
RectType: TYPE = ObjectCast.RectType;-- {up, down, front, back, left, right};
TubeSurface: TYPE = REF TubeSurfaceObj;
TubeSurfaceObj: TYPE = ObjectCast.TubeSurfaceObj;
DiskSurface: TYPE = REF DiskSurfaceObj;
DiskSurfaceObj: TYPE = ObjectCast.DiskSurfaceObj;
DiskType: TYPE = ObjectCast.DiskType;-- {top, bottom};
ShellSurface: TYPE = REF ShellSurfaceObj;
ShellSurfaceObj: TYPE = ObjectCast.ShellSurfaceObj;
Composite: TYPE = REF CompositeObj;
CompositeObj: TYPE = CSG.CompositeObj;
CSGTree: TYPE = REF CSGTreeObj;
CSGTreeObj: TYPE = CSG.CSGTreeObj;
PointSetOp: TYPE = CSG.PointSetOp;
Classification: TYPE = REF ClassificationObj;
ClassificationObj: TYPE = CastRays.ClassificationObj;
Ray: TYPE = REF RayObj;
RayObj: TYPE = CastRays.RayObj;
HitRec: TYPE = REF HitRecObj;
HitRecObj: TYPE = ObjectCast.HitRecObj;
HitArray: TYPE = REF HitArrayObj;
HitArrayObj: TYPE = ObjectCast.HitArrayObj;
SurfaceArray: TYPE = REF SurfaceArrayObj;
SurfaceArrayObj: TYPE = CSG.SurfaceArrayObj;
ParameterArray: TYPE = CastRays.ParameterArray;
InOutArray: TYPE = CastRays.InOutArray;
NormalArray: TYPE = CastRays.NormalArray;
globalHitArray: HitArray;
globalHitArray2: HitArray;
BlockCast:
PUBLIC
PROC [localRay: Ray, surfs: SurfaceArray, prim: Primitive]
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;
hitCount: NAT ← 0;
doesHit: BOOL;
upS, downS, frontS, backS, rightS, leftS: Surface;
almostZero: REAL ← 1.0e-12;
upS ← surfs[1]; downS ← surfs[2]; frontS ← surfs[3];
backS ← surfs[4]; rightS ← surfs[5]; leftS ← surfs[6];
class ← CastRays.GetClassFromPool[];
p1d2 ← localRay.basePt[1]*localRay.direction[2];
p1d3 ← localRay.basePt[1]*localRay.direction[3];
p2d1 ← localRay.basePt[2]*localRay.direction[1];
p2d3 ← localRay.basePt[2]*localRay.direction[3];
p3d1 ← localRay.basePt[3]*localRay.direction[1];
p3d2 ← localRay.basePt[3]*localRay.direction[2];
front face and back face
IF Abs[localRay.direction[3]] > almostZero
THEN {
D ← Abs[localRay.direction[3]];
front face
xd ← p1d3 + localRay.direction[1] - p3d1;
IF -
D <= xd
AND xd <=
D
THEN {
yd ← p2d3 + localRay.direction[2] - p3d2;
IF -D <= yd AND yd <= D
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit
THEN {
t ← (1.0 - localRay.basePt[3])/localRay.direction[3];
IF t>0
THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, frontS, [0,0,1]];
};
};
back face
xd ← p1d3 - localRay.direction[1] - p3d1;
IF -
D <= xd
AND xd <=
D
THEN {
yd ← p2d3 - localRay.direction[2] - p3d2;
IF -D <= yd AND yd <= D
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit
THEN {
t ← (-1.0 - localRay.basePt[3])/localRay.direction[3];
IF t>0
THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, backS, [0,0,-1]];
};
};
};
up face and down face
IF Abs[localRay.direction[2]] > almostZero
THEN {
D ← Abs[localRay.direction[2]];
up face
xd ← p1d2 + localRay.direction[1] - p2d1;
IF -
D <= xd
AND xd <=
D
THEN {
zd ← p3d2 + localRay.direction[3] - p2d3;
IF -D <= zd AND zd <= D
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit
THEN {
hitCount ← hitCount + 1;
t ← (1.0 - localRay.basePt[2])/localRay.direction[2];
hits[hitCount]^ ← [t, upS, [0,1,0]]};
down face
xd ← p1d2 - localRay.direction[1] - p2d1;
IF -
D <= xd
AND xd <=
D
THEN {
zd ← p3d2 - localRay.direction[3] - p2d3;
IF -D <= zd AND zd <= D
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit
THEN {
t ← (-1.0 - localRay.basePt[2])/localRay.direction[2];
IF t>0
THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, downS, [0,-1,0]];
};
};
};
left face and right face
IF Abs[localRay.direction[1]] > almostZero
THEN {
D ← Abs[localRay.direction[1]];
right face
yd ← p2d1 + localRay.direction[2] - p1d2;
IF -
D <= yd
AND yd <=
D
THEN {
zd ← p3d1 + localRay.direction[3] - p1d3;
IF -D <= zd AND zd <= D
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit
THEN {
t ← (1.0 - localRay.basePt[1])/localRay.direction[1];
IF t>0
THEN {
hitCount ← hitCount + 1;
hits[hitCount]^ ← [t, rightS, [1,0,0]];
};
};
left face
yd ← p2d1 - localRay.direction[2] - p1d2;
IF -
D <= yd
AND yd <=
D
THEN {
zd ← p3d1 - localRay.direction[3] - p1d3;
IF -D <= zd AND zd <= D
THEN doesHit ← TRUE ELSE doesHit ← FALSE}
ELSE doesHit ← FALSE;
IF doesHit
THEN {
t ← (-1.0 - localRay.basePt[1])/localRay.direction[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
Count as a miss (for now).
class.count ← 0;
class.classifs[1] ← FALSE;
class.count ← 1;
class.params[1] ← hits[1].t;
class.surfaces[1] ← hits[1].surf;
class.classifs[1] ← FALSE; class.classifs[2] ← FALSE;
class.normals[1] ← hits[1].normal;
globalOutHandle.PutChar['#];
};
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;
globalOutHandle.PutChar['*];
};
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 [rootArray: ARRAY[1..2] OF REAL, rootCount: NAT] RETURNS [posRootArray: ARRAY[1..2] OF REAL, posRootCount: NAT] = {
Use only the positive roots.
i: NAT ← 1;
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;
};
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 = x0 + t*Dx
y = y0 + t*Dy
z = z0 + t*Dz
hence, at the intersection of ray and sphere, t satisfies:
(x0 + t*Dx)^2 + (y0 + t*Dy)^2 + (z0 + t*Dz)^2 = r^2
which reduces to:
t^2 * (Dx^2 + Dy^2 + Dz^2) + t * 2(x0*Dx + y0*Dy + z0*Dz) +
(x0^2 + y0^2 + z0^2 - r^2) = 0;
Let p = [x0, y0, z0], and d = [Dx, Dy, Dz] then in dot product notation:
t^2 * d.d + t * 2 * d.p + p.p - r^2 = 0. Use quadratic formula
d: Vector ← localRay.direction;
p: Point3d ← localRay.basePt;
tArray: ARRAY[1..2] OF REAL;
rootCount: NAT;
dp: REAL ← SVVector3d.DotProduct[d,p];
dp2: REAL ← dp*dp;
dd: REAL ← SVVector3d.DotProduct[d,d];
pp: REAL ← SVVector3d.DotProduct[p,p];
t1, t2: REAL;
class ← CastRays.GetClassFromPool[];
[tArray, rootCount, ----] ← PositiveRoots[dd,2*dp,pp-r*r];
SELECT rootCount
FROM
0 => {
class.count ← 0;
class.classifs[1] ← FALSE;
};
1 => {
Treat as a miss for now.
class.count ← 0;
class.classifs[1] ← FALSE;
class.count ← 1;
class.classifs[1] ← FALSE;
class.classifs[2] ← FALSE;
class.params[1] ← tArray[1];
class.surfaces[1] ← thisShell;
class.normals[1][1] ← p[1] + d[1]*class.params[1];
class.normals[1][2] ← p[2] + d[2]*class.params[1];
class.normals[1][3] ← p[3] + d[3]*class.params[1];
class.primitives[1] ← prim;
};
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]*class.params[1];
class.normals[1][2] ← p[2] + d[2]*class.params[1];
class.normals[1][3] ← p[3] + d[3]*class.params[1];
class.normals[2][1] ← p[1] + d[1]*class.params[2];
class.normals[2][2] ← p[2] + d[2]*class.params[2];
class.normals[2][3] ← p[3] + d[3]*class.params[2];
};
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)
d: Vector ← localRay.direction;
p: Point3d ← localRay.basePt;
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: Vector;
p1d2, p2d1, p2d3, p3d2, d2d2, xd2, zd2: REAL;
topS ← surfaces[1]; bottomS ← surfaces[2]; tubeS ← surfaces[3];
class ← CastRays.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 (or perhaps this is an error). Treat as a miss for now.
class.count ← 0;
class.classifs[1] ← FALSE;
};
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.
class.count ← 0;
class.classifs[1] ← FALSE;
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 shoud 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, count the ray as a miss:
IF tArray[1] <= 0
THEN {
class.count ← 0;
class.classifs[1] ← FALSE;
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. One on a disk face.
SELECT planeHitCount FROM
0 => {
-- may be skimming a rim and only the tube notices
Count as a miss for now
class.count ← 0;
class.classifs[1] ← FALSE;
};
1 => {
If either the tubeHit or the planeHit is negative, the ray originates from inside of the cylinder. Count as a miss for now.
IF tubeHits[1].t <= 0
OR planeHits[1].t <=0
THEN {
class.count ← 0;
class.classifs[1] ← FALSE;
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. Count as a miss for now.
IF planeHits[1].t <=0 OR planeHits[2].t <= 0 THEN {
class.count ← 0;
class.classifs[1] ← FALSE;
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
class.count ← 0;
class.classifs[1] ← FALSE;
};
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.
Treat as a miss for now.
class.count ← 0;
class.classifs[1] ← FALSE;
};
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 OR planeHits[2].t <= 0 THEN {
class.count ← 0;
class.classifs[1] ← FALSE;
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..ObjectCast.globalObjDepth]
DO
globalHitArray[i] ← NEW[HitRecObj];
globalHitArray2[i] ← NEW[HitRecObj];
ENDLOOP;
};
Init[];