DIRECTORY
Algebra3d,
CastRays,
ObjectCast,
CSG,
DisplayList3d,
Matrix3d,
Rope,
SVVector3d;

ObjectCastImpl: PROGRAM
IMPORTS Algebra3d, CastRays, SVVector3d
EXPORTS ObjectCast =

BEGIN

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] = {

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];

IF Abs[localRay.direction[3]] > almostZero THEN {
D _ Abs[localRay.direction[3]];

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]];
};
};

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]];
};
};
};

IF Abs[localRay.direction[2]] > almostZero THEN {
D _ Abs[localRay.direction[2]];

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]]};

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]];
};
};
};

IF Abs[localRay.direction[1]] > almostZero THEN {
D _ Abs[localRay.direction[1]];

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]];
};
};

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
class.count _ 0;
class.classifs[1] _ FALSE;
};
2 => { -- ray enters and exits block properly
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;
};
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] = {
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] = {
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] = {
r: REAL = 1.0;
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 => {
class.count _ 0;
class.classifs[1] _ FALSE;
};
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] = {
yLow: REAL = -1.0;  yHigh: REAL = 1.0;
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;
tubeNormal, tempNormal: Vector;
p1d2, p2d1, p2d3, p3d2, d2d2, xd2, zd2: REAL;

topS _ surfaces[1];  bottomS _ surfaces[2];  tubeS _ surfaces[3];
class _ CastRays.GetClassFromPool[];

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];

IF Abs[d[2]] > almostZero THEN {

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] ];
};
};

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
};

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=> {
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 => {
class.count _ 0;
class.classifs[1] _ FALSE;
};
2 => {
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.
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
class.count _ 0;
class.classifs[1] _ FALSE;
};
1 => {
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 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.
class.count _ 0;
class.classifs[1] _ FALSE;
};
2 => { -- goes through both disk faces
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
};
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 = {
globalHitArray _ NEW[HitArrayObj];
globalHitArray2 _ NEW[HitArrayObj];
FOR i: NAT IN[1..ObjectCast.globalObjDepth] DO
globalHitArray[i] _ NEW[HitRecObj];
globalHitArray2[i] _ NEW[HitRecObj];
ENDLOOP;
};

Init[];

END.



¦File: ObjectCastImpl.mesa
Author: Eric Bier on July 20, 1982 3:36 pm
Last edited by Bier on August 4, 1983 5:24 pm
Contents: Ray-Object intersection tests for block, sphere, and cylinder
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
front face and back face
front face
back face
up face and down face
up face
down face
left face and right face
right face
left face
Count as a miss (for now).
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['#];
The hit with the lower value of t is the entry point.  The other is the exit.
globalOutHandle.PutChar['*];
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;
};
Use only the positive roots.  Preserves the order of the roots, so they will still be in increasing order.
Count the positive roots.  Preserves the order of the roots, so they will still be in increasing order.
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.
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

Treat as a miss for now.
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 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

r: REAL = 1.0;(built into code)
thisShell: Surface;
Do tube test.
Do 2 ray-plane tests

.top face and bottom face
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.
bottom face.
Sort the plane hits.
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.
Ray just scratches tube.  Count as a miss.
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;
};
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.
At this point, if any of the intersections are on the wrong side of the ray, count the ray as a miss:
Count as a miss for now
If either the tubeHit or the planeHit is negative, the ray originates from inside of the cylinder.  Count as a miss for now.
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.
Treat as a miss for now.
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.
End of zero tube hits of 2 tube roots.
Create a global HitArray.
�Ê¡–

"cedar" style˜Iheadš
œ™Iprocš
œ*™*Lš
œ-™-Lš
œG™GL˜š
Ïk	˜	Lš
œ
˜
Lš
œ	˜	Lš
œ˜Lš�œ
˜Lš
œ˜Lš
œ	˜	Lš
œ˜Lš
œ˜—L˜šœ�˜Lš�œ ˜'Lš�œ
˜—L˜Lš
�˜˜Lšœ	�œ˜!Lšœ	�œ˜!Lšœ�œ˜'Lšœ�œ˜!Lšœ	�œ�œ
�œ
˜Lšœ�œ!˜6L˜Lšœ�œ�œ˜#Lšœ�œ�œ˜&L˜Lšœ
�œ�œ˜'Lšœ�œ˜1Lšœ
�œÏc(˜MLšœ
�œ�œ˜'Lšœ�œ˜1Lšœ
�œ�œ˜'Lšœ�œ˜1Lšœ
�œž˜6Lšœ�œ�œ˜)Lšœ�œ˜3L˜Lšœ�œ�œ˜#Lšœ�œ�œ˜&L˜Lšœ	�œ�œ˜Lšœ�œ�œ˜"L˜Lšœ�œ�œ˜"L˜Lšœ�œ�œ˜.Lšœ�œ˜5L˜Lšœ�œ�œ˜Lšœ�œ˜L˜Lšœ�œ�œ˜Lšœ�œ˜'Lšœ
�œ�œ
˜!Lšœ
�œ˜+L˜Lšœ�œ�œ˜)Lšœ�œ�œ˜,L˜Lšœ�œ˜/Lšœ�œ˜'Lšœ
�œ˜)L˜Lš
œ˜Lš
œ˜L˜L˜šÏn	œ�œ
�œ7�œ˜pLš
œÿ
™ÿ
Lš
œ¸™¸Lš
œ�™�Lš
œŠ™ŠLš
œ¡™¡L˜Lšœ�œ
˜Lšœ�
œ�œ
˜Lš
œ ˜ Lšœ$�œ
˜)Lšœ
�œ˜Lšœ	�œ
˜Lš
œ2˜2Lšœ�œ˜Lš
œ6˜6Lš
œ8˜8Lš
œ$˜$Lš
œ0˜0Lš
œ0˜0Lš
œ0˜0Lš
œ0˜0Lš
œ0˜0Lš
œ0˜0L˜Lš
œ™š�œ)�œ˜1Lš�
œ˜ Lš
œ
™
Lš
œ)˜)š
�œ�
œ�œ�
œ
�œ˜Lš
œ)˜)Lš�œ�
œ�œ�
˜Lš�œ�œ
�œ�œ
˜)—Lš�œ�œ
˜š�œ	�œ˜Lš
œ5˜5š�œ�œ˜
Lš
œ˜Lš
œ'˜'Lš
œ˜—Lš
œ˜—Lš
œ	™	Lš
œ)˜)š
�œ�
œ�œ�
œ
�œ˜Lš
œ)˜)Lš�œ�
œ�œ�
˜Lš�œ�œ
�œ�œ
˜)—Lš�œ�œ
˜š�œ	�œ˜Lš
œ6˜6š�œ�œ˜
Lš
œ˜Lš
œ'˜'Lš
œ˜—Lš
œ˜—Lš
œ˜—L˜Lš
œ™š�œ)�œ˜1Lš�
œ˜ Lš
œ™Lš
œ)˜)š
�œ�
œ�œ�
œ
�œ˜Lš
œ)˜)Lš�œ�
œ�œ�
˜Lš�œ�œ
�œ�œ
˜)—Lš�œ�œ
˜š�œ	�œ˜Lš
œ˜Lš
œ5˜5Lš
œ&˜&—Lš
œ	™	Lš
œ)˜)š
�œ�
œ�œ�
œ
�œ˜Lš
œ)˜)Lš�œ�
œ�œ�
˜Lš�œ�œ
�œ�œ
˜)—Lš�œ�œ
˜š�œ	�œ˜Lš
œ6˜6š�œ�œ˜
Lš
œ˜Lš
œ'˜'Lš
œ˜—Lš
œ˜—Lš
œ˜—L˜Lš
œ™š�œ)�œ˜1Lš�
œ˜ Lš
œ
™
Lš
œ)˜)š
�œ�
œ�œ�
œ
�œ˜Lš
œ)˜)Lš�œ�
œ�œ�
˜Lš�œ�œ
�œ�œ
˜)—Lš�œ�œ
˜š�œ	�œ˜Lš
œ5˜5š�œ�œ˜
Lš
œ˜Lš
œ'˜'Lš
œ˜—Lš
œ˜—Lš
œ	™	Lš
œ)˜)š
�œ�
œ�œ�
œ
�œ˜Lš
œ)˜)Lš�œ�
œ�œ�
˜Lš�œ�œ
�œ�œ
˜)—Lš�œ�œ
˜š�œ	�œ˜Lš
œ6˜6š�œ�œ˜
Lš
œ˜Lš
œ'˜'Lš
œ˜—Lš
œ˜—Lš
œ˜—L˜š�œ
�˜šœž˜&Lš
œ˜Lšœ�œ
˜Lš
œ˜—šœžG˜NLš
œ™Lš
œ˜Lšœ�œ
˜Lš
œ™Lš
œ™Lš
œ!™!Lšœ�œ�œ
™6Lš
œ"™"Lš
œ™Lš
œ˜—šœž&˜-Lš
œM™MLš
œ˜š�œ�œ˜Lš
œ˜Lš
œ˜Lš
œ!˜!Lš
œ!˜!Lš
œ"˜"Lš
œ#˜#—š�œ˜Lš
œ˜Lš
œ˜Lš
œ!˜!Lš
œ!˜!Lš
œ"˜"Lš
œ$˜$—Lšœ�œ�œ�œ
˜QLš
œ4˜4Lš
œ˜—šœžŠ
˜‘
Lš
œ<˜<Lš
œ˜Lš
œ"˜"Lš
œ"˜"Lš
œ'˜'Lš
œ'˜'Lš
œ(˜(Lš
œ(˜(Lšœ�œ�œ�œ˜RLš
œ2˜2Lš
œ™Lš
œ˜—Lš�œ�œ˜ —Lšœž˜—L˜Lšœ�œ�œ
˜Lš	œ
�œ�œ�œ
�œ˜)šŸœ�œ�œ�œ
�œ�œ˜CLšœ�œ
˜Lšœ�œ
˜Lš
œ˜š�œ�œ
�œ�˜š�œ�œ
�œ	�˜š�œ�œ˜Lš
œ
˜
Lš
œ˜Lš
œ˜Lš
œ˜Lš
œ˜Lš
œ˜——Lš�œ
˜—Lš�œ
˜Lšœž˜J˜—šŸ
œ�	œ
�œ�œ
�œ
�œ�œ�œ�œ
�œ�œ™Œ
J™Lšœ�œ™Lš�œ�™š�œ�œ™š�œ�œ
�œ�™!Lš
œ™—Lš�œ
™Lš
œ™L™
—Lš�œ
™Lš�œ
™Jš
œ™—šœ
Ÿ
œ�	œ�œ�œ
�œ�œ
�œ�œ˜vJ™jLšœ�œ˜Lš
œ=˜=Lš�œ�œ
�œ
˜š�œ�˜š�œ�œ˜š�œ�œ
�œ�˜!Lš
œ˜—Lš�œ
˜Lš
œ˜L˜
—Lš�œ
˜—Lš�œ
˜Jš
œ˜—šŸœ�	œ�œ�œ
�œ�œ
�œ%�œ˜xJ™gLš
œ˜Lš
œB˜BLš�œ�œ
�œ
˜"š�œ�œ
�œ�˜$Lš�œ�œ+˜D—Lš�œ
˜Jš
œ˜—š	œ
Ÿ
œ�œ
�œ6�œ˜qLš
œ™Lš
œ9™9Lš
œ4™4Lšœ�œ˜Lš
œ$™$Lš
œ™Lš
œ™Lš
œ:™:Lš
œ3™3Lš
œ™Lš
œ;™;Lš
œ™Lš
œH™HLš
œ?™?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š
œ2™2Lš
œ2™2Lš
œ2™2Lš
œ™Lš
œ˜—š
œ˜Lš
œ˜Lš
œM˜MLš
œ,˜,Lš
œ=˜=Lšœ�œ�œ�œ
˜QLš
œ1˜1Lš
œ2˜2Lš
œ2˜2Lš
œ2˜2Lš
œ2˜2Lš
œ2˜2Lš
œ2˜2Lš
œ˜—Lš�œ�œ
˜—Lš
œ˜—L˜šŸœ�œ
�œ:�œ˜vLš
œ%™%Lš
œ8™8Lš
œ'™'Lš
œ™Lš
œ™Lš
œ™Lš
œ9™9Lš
œ™Lš
œJ™JLš
œQ™QLš
œ™Lšœ�œ�œ˜&Lšœ�œ™Lš
œ˜Lš
œ˜Lšœ�œ�œ
�œ
˜Lšœ	�œ
ž"˜0Lšœ#�œ
˜'Lšœ�œ˜L˜Lšœ�œ
˜Lš
œ$˜$Lšœ�œ˜Lš
œ&˜&Lšœ�œ˜Lšœ�œ
˜Lš
œ+˜+Lš
œ™Lš
œ˜Lšœ(�œ
˜-L˜Lš
œA˜ALš
œ$˜$L˜—š
œ
™
Lš
œ˜Lš
œ˜Lšœ ž˜,Lš
œA˜A—Lš
œ™š
œ™Lš�œ�œ˜!š
œ	™	Lš
œ>™>Lš
œ7™7Lš
œU™ULš
œY™YLš
œO™OLš
œR™RLš
œ1™1Lš
œ!™!Lš
œ¯
™¯
Lš
œ˜Lš
œ˜Lš
œ˜Lš
œ˜Lš
œ˜Lš
œ˜Lš
œ˜Lš
œ˜Lš�œž
˜!š�œ˜Lš
œ˜š�œ�œ˜Lš
œ"˜"Lš
œ2˜2Lš
œ˜—Lš
œ˜——š
œ™Lš
œ˜Lš
œ˜Lš
œ˜Lš�œž
˜!š�œ˜Lš
œ˜š�œ�œ˜Lš
œ"˜"Lš
œ5˜5Lš
œ˜—Lšœž$˜'—Lš
œ˜—L™š�œ�œ˜š�œ!�œ˜)Lš
œ˜Lš
œ ˜ Lš
œ!˜!Lš
œ ˜ Lš
œ&˜&Lš
œ*˜*Lš
œ˜Lš
œ ˜ Lš
œ!˜!Lš
œ˜—Lš
œ˜—Lš
œ.˜.š�œ�˜šœžD˜Jš�œ�˜š
œ˜Lš
œ˜Lšœ�œ
˜Lš
œ˜—š
œ˜Lšœ
ž�
™‚
Lš
œ˜Lšœ�œ
˜Lš
œ˜—š
œ˜Lš
œ˜Lš
œL˜LLš
œ!˜!Lš
œ&˜&Lš
œ&˜&Lš
œ'˜'Lš
œ'˜'Lšœ�œ�œ�œ
˜QLš
œ5˜5Lšœž˜—Lš�œ�œ
˜—Lšœž˜—š
œ˜Lšœ
ž)™*Lš
œ˜Lšœ�œ
˜Lš
œ™š�œ�œ�œž™>Lš�œ�œ
�œž*™KLš
œ™Lšœ�œ
™Lšœ�œ
™Lš
œ™Lš
œ™Lš
œ™Lš
œK™KLš
œ,™,Lš
œ™Lš
œ,™,Lš
œ
™
—š�œž	™Lšœ�œ
™Lš
œ™—Lš
œ˜—Lš
œ˜Lš
žp™pLš
œ1™1Lš
œ(™(Lš
œT™Tš
œ7™7Lš
œ˜Lš
œ˜š�œ�œ�œž˜=Lš
œ ˜ Lš
œ&˜&Lš
œ˜Lš
œ&˜&Lš
œ:˜:—Lš
œ˜š�œ�œ�œž˜>Lš
œ ˜ Lš
œ&˜&Lš
œ˜Lš
œ&˜&Lš
œ:˜:—š�œ�˜šœž%˜,Lš�œ�œ
�œž|˜œ
L™eš�œ�œ˜Lš
œ˜Lšœ�œ
˜Lš�œ
˜Lš
œ˜—Lš
œ˜Lš
œN˜NLš
œ ˜ Lš
œ&˜&Lš
œ&˜&Lš
œ5˜5Lšœ�œ�œ�œ
˜QLš
œ5˜5Lš
œ˜—šœž:˜ALš�œ�˜šœž2˜9Lš
œ™Lš
œ˜Lšœ�œ
˜Lš
œ˜—š
œ˜L™|š�œ�œ�œ˜2Lš
œ˜Lšœ�œ
˜Lš�œ
˜Lš
œ˜—Lš
œ˜š�œ �œ˜(Lš
œ ˜ Lš
œ!˜!Lš
œ%˜%Lš
œ&˜&Lš
œ&˜&Lš
œ(˜(—š�œ˜Lš
œ ˜ Lš
œ!˜!Lš
œ%˜%Lš
œ&˜&Lš
œ&˜&Lš
œ)˜)—Lšœ�œ�œ�œ
˜QLš
œ5˜5Lš
œ˜—šœž8˜>L™†
š
œ3˜3Lš
œ˜Lšœ�œ
˜Lš�œ
˜Lš
œ˜—Lš
œ˜š�œ!�œ˜)Lš
œ!˜!Lš
œ!˜!Lš
œ&˜&Lš
œ&˜&Lš
œ'˜'Lš
œ(˜(—š�œ˜Lš
œ!˜!Lš
œ!˜!Lš
œ&˜&Lš
œ&˜&Lš
œ'˜'Lš
œ)˜)—Lšœ�œ�œ�œ
˜QLš
œ5˜5Lš
œ˜—Lš�œ�œ
ž'˜8Lš
œ˜—šœž,˜2š�œ�˜šœž	˜Lš
œ˜Lšœ�œ
˜Lš
œ˜—šœž•
˜œ
Lš
œ™Lš
œ˜Lšœ�œ
˜Lš
œ˜—šœž˜&L™†
š
œ3˜3Lš
œ˜Lšœ�œ
˜Lš�œ
˜Lš
œ˜—Lš
œ˜š�œ!�œ˜)Lš
œ!˜!Lš
œ!˜!Lš
œ&˜&Lš
œ&˜&Lš
œ'˜'Lš
œ(˜(—š�œ˜Lš
œ!˜!Lš
œ!˜!Lš
œ&˜&Lš
œ&˜&Lš
œ'˜'Lš
œ)˜)—Lšœ�œ�œ�œ
˜QLš
œ5˜5Lš
œ˜—Lš�œ�œž$˜6Lš
œ&™&L˜——Lš�œ�œžœž˜D—Lšœž˜—Lš�œ�œ
ž)˜:—Lšœž!˜#—L˜š
Ÿœ�œ�œ�œ�œ˜&Lš�œ�œ�œ�œ˜ Lš
œ˜—L˜šŸœ�œ˜Lš
œ™Lšœ�œ˜"Lšœ�œ˜#š�œ�œ
�œ�˜.Lšœ�œ˜#Lšœ�œ˜$—Lš�œ
˜Lš
œ˜—L˜Lš
œ˜L˜—Lš�œ
˜L˜L˜L˜—…—@Æt