File: ObjectCastImplC.mesa
Last edited by Bier on January 2, 1983 8:31 pm
Author: Bier on July 4, 1983 10:44 pm
DIRECTORY
BasicObject3d,
CastRays,
CSG,
ObjectCast,
Polynomial,
RealFns,
SV2d,
Roots;
ObjectCastImplC: PROGRAM
IMPORTS CastRays, Polynomial, RealFns
EXPORTS ObjectCast =
BEGIN
A torus has a single surface. Classifying a ray with respect to this surface involves solving a single quadric equation.
Let the radius from the center of the doughnut hole to the center of the doughnut be C.
Let the cross-sectional radius be B.
Let the doughnut be a circle revolved around the y axis.
Consider a cross section of the doughnut. Call the line parallel to the y axis and running through the center of the cross section L. Let r(y) be the distance perpendicular to L from L to the circular boundary of the cross section. Then:
(1) r^2 + y^2 = C^2
From which we find:
(2) r = +- SqRt[C^2 - y^2]
Thus the two concentric circles cut by the plane y = y0 have equations:
x^2 + z^2 = (B +- SqRt[C^2 - y^2])^2 or compactly:
(3) x^2 + z^2 = (B +- r)^2 = B^2 +- 2Br + r^2
We can get rid of the square root by regrouping:
(4) x^2 + z^2 - B^2 - r^2 = +- 2Br
and squaring:
(5) (x^2 + z^2 - B^2 - r^2)^2 = 4B^2r^2
Expanding, we have:
x^4 + x^2z^2 - x^2*B^2 - x^2*r^2 +
z^2*x^2 + z^4 - z^2*B^2 - z^2*r^2
-B^2*x^2 -B^2*z^2 + B^4 + B^2*r^2
-r^2*x^2 - r^2*z^2 + B^2*r^2 + r^4 = 4*B^2*r^2;  OR

x^4 + z^4 + B^4 + r^4 +
2*x^2z^2 - 2*x^2*B^2 - 2*x^2*r^2 - 2*z^2*B^2 - 2*z^2*r^2 + 2*B^2*r^2
 = 4*B^2*r^2

Plugging in r^2 = (C^2 - y^2), r^4 = (C^4 - 2*C^2*y^2 + y^4):

x^4 + z^4 + B^4 + C^4 - 2*C^2*y^2 + y^4 +
2*x^2z^2 - 2*x^2*B^2 - 2*x^2*(C^2 - y^2) - 2*z^2*B^2 - 2*z^2*(C^2 - y^2) + 2*B^2*(C^2 - y^2)
 = 4*B^2*(C^2 - y^2) = 4*B^2*C^2 - 4*B^2*y^2 OR

x^4 + y^4 + z^4 + B^4 + C^4 +
2*x^2z^2 + 2*x^2*y^2 + 2*z^2*y^2
x^2 (- 2*B^2 - 2*C^2) +
y^2 (- 2*B^2- 2*C^2) +
z^2 (- 2*B^2 - 2*C^2) +
- 2*C^2 + 2*B^2*C^2 = 4*B^2*C^2 - 4*B^2*y^2 OR

x^4 + y^4 + z^4 + B^4 + C^4 +
2*x^2z^2 + 2*x^2*y^2 + 2*z^2*y^2
x^2 (- 2*B^2 - 2*C^2) +
y^2 (2*B^2 - 2*C^2) +
z^2 (- 2*B^2 - 2*C^2) +
- 2*B^2*C^2 = 0

Let D = B^2 + C^2, E = B^2 - C^2, F = C^2B^2 Then:

(6) x^4 + y^4 + z^4 + B^4 + C^4 +
  + 2(x^2*y^2 + x^2z^2 + y^2z^2
   - Dx^2 + Ey^2 - Dz^2 - F) = 0


Now recall the ray equations (using single letters for all variables for compactness):
(7) x = n + t*N
(8) y = q + t*Q
(9) z = l + t*L
We wish to solve for t by plugging (7) - (9) into (6). We get:

(n + t*N)^4 + (q + t*Q)^4 + (l + t*L)^4 + B^4 + C^4
+ 2( (n + t*N)^2*(q + t*Q)^2 + (n + t*N)^2*(l + t*L)^2 + (q + t*Q)^2*(l + t*L)^2
 -D(n + t*N)^2 + E(q + t*Q)^2 - D(l + t*L)^2 - F) = 0

If I've done my algebra correctly, this boils down to:
a4*t^4 + a3*t^3 + a2*t^2 + a1*t + a0
Where:
a4 = N^4 + Q^4 + L^4 + 2(N^2(Q^2+L^2) + Q^2L^2);
a3 = 4(nN^3 + qQ^3 + lL^3 + nN(Q^2 + L^2) + N^2(qQ + lL) + lLQ^2 + qQL^2)
a2 = 2(3n^2N^2 + 3q^2Q^2 + 3l^2L^2 + N^2(q^2 + l^2) + n^2(Q^2+L^2) + 4nN(qQ + lL)
- DN^2 + EQ^2 - DL^2 + l^2Q^2 + L^2q^2 + 4qQlL)
a1 = 4(n^3N + q^3Q + l^3L + nN(q^2 + l^2) + (qQ + lL)n^2 - nND + qQE - lLD + qQl^2 + lLq^2)
a0 = n^4 + q^4 + l^4 + B^4 + C^4 +
2(n^2(q^2 + l^2) - Dn^2 + Eq^2 - Dl^2+ q^2l^2 - F)
Here goes:
Classification: TYPE = CSG.Classification;
Point2d: TYPE = SV2d.Point2d;
Primitive: TYPE = CSG.Primitive;
Ray: TYPE = CSG.Ray;
Surface: TYPE = CSG.Surface;
TorusRec: TYPE = BasicObject3d.TorusRec;
globalPoly: Polynomial.Ref;
PositiveRoots: PROCEDURE [a4, a3, a2, a1, a0: REAL] RETURNS [rootArray: REF Polynomial.RealRootRec] = {
Use only the positive roots.
i: NAT ← 1;
globalPoly[4] ← a4; -- optional
globalPoly[3] ← a3; -- optional
globalPoly[2] ← a2; -- optional
globalPoly[1] ← a1; -- optional
globalPoly[0] ← a0; -- optional
rootArray ← Polynomial.RealRoots[globalPoly]; -- optional
[realRootsOneOrigin, rootCount] ← Roots.RealQuartic[a4, a3, a2, a1, a0];
realRoots[0] ← realRootsOneOrigin[1]; realRoots[1] ← realRootsOneOrigin[2];
realRoots[2] ← realRootsOneOrigin[3]; realRoots[3] ← realRootsOneOrigin[4];
IF rootArray.nRoots = 0 THEN RETURN;
WHILE i <= rootArray.nRoots DO
IF rootArray[i-1] <= 0 THEN {
FOR j: NAT IN [i..rootArray.nRoots-1] DO
rootArray[j-1] ← rootArray[j];
ENDLOOP;
rootArray.nRoots ← rootArray.nRoots - 1;
}
ELSE {i ← i + 1};
ENDLOOP;
};

ToroidCast
: PUBLIC PROC [localRay: Ray, prim: Primitive, torusRec: REF ANY, surface: Surface]
RETURNS [class: Classification] = {
torusData: TorusRec;
realRoots: REF Polynomial.RealRootRec; -- optional
realRoots: ARRAY[0..3] OF REAL;
realRootsOneOrigin: ARRAY[1..4] OF REAL;
rootCount: NAT;
B,C: REAL;
n,N,q,Q,l,L: REAL;
BB, CC: REAL;
D,E,F: REAL;
nn, NN, qq, QQ, ll, LL, NNNN, QQQQ, LLLL, nN, qQ, lL: REAL;
a4,a3,a2,a1,a0: REAL;
class ← CastRays.GetClassFromPool[];
torusData ← NARROW[torusRec];
B ← torusData.bigRadius;
C ← torusData.sectionRadius;
n ← localRay.basePt[1]; N ← localRay.direction[1];
q ← localRay.basePt[2]; Q ← localRay.direction[2];
l ← localRay.basePt[3]; L ← localRay.direction[3];
BBB*B;
CCC*C;
DBB + CC;
EBB - CC;
FBB*CC;
nn ← n*n; NNN*N; qq ← q*q; QQQ*Q; ll ← l*l; LLL*L;
NNNNNN*NN; QQQQQQ*QQ; LLLLLL*LL;
nN ← n*N; qQ ← q*Q; lL ← l*L;
a4 ← NNNN + QQQQ + LLLL + 2*(NN*(QQ + LL) + QQ*LL);
a3 ← 4*(nN*NN + qQ*QQ + lL*LL + nN*(QQ + LL) + NN*(qQ + lL) + lL*QQ + qQ*LL);
a2 ← 2*(3*(nn*NN + qq*QQ + ll*LL) + NN*(qq + ll) + nn*(QQ+LL) + 4*nN*(qQ+lL)
  - D*NN + E*QQ - D*LL + ll*QQ + LL*qq + 4*qQ*lL);
a1 ← 4*(nn*nN + qq*qQ + ll*lL + nN*(qq + ll) + (qQ + lL)*nn - nN*D + qQ*E - lL*D + qQ*ll + lL*qq);
a0 ← nn*nn + qq*qq + ll*ll + BB*BB + CC*CC +
 2*(nn*(qq + ll) - D*nn + E*qq - D*ll + qq*ll - F);
realRoots ← PositiveRoots[a4, a3, a2, a1, a0];
SELECT realRoots.nRoots FROM
SELECT rootCount FROM
0 => {-- complete miss.
class.count ← 0;
class.classifs[1] ← FALSE;
};
1 => {-- a skim. Count as a miss.
class.count ← 0;
class.classifs[1] ← FALSE;
};
2 => {-- cuts one cross section of the doughnut. Sort roots by size.
x, y, z, theta: REAL;
class.count ← 2;
IF realRoots[0] < realRoots[1] THEN
{class.params[1] ← realRoots[0];
class.params[2] ← realRoots[1]}
ELSE
{class.params[2] ← realRoots[0];
class.params[1] ← realRoots[1]};
class.surfaces[1] ← surface; class.surfaces[2] ← surface;
class.classifs[1] ← FALSE; class.classifs[2] ← TRUE; class.classifs[3] ← FALSE;
class.primitives[1] ← class.primitives[2] ← prim;
Calculate the surface normals at the two hits points.
x ← n + class.params[1]*N; y ← q + class.params[1]*Q; z ← l + class.params[1]*L;
theta ← RealFns.ArcTan[-z,x];
class.normals[1] ← [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z];
To normalize, divide by C. (I won't do this since the lighting model normalizes anyway)
x ← n + class.params[2]*N; y ← q + class.params[2]*Q; z ← l + class.params[2]*L;
theta ← RealFns.ArcTan[-z,x];
class.normals[2] ← [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z];
};
3 => {-- cuts one cross section. Skims the other. For now, count as a miss.
class.count ← 0;
class.classifs[1] ← FALSE;
};
4 => {-- cuts both cross sections. Sort the roots. Pair them in order.
x, y, z, theta: REAL;
SortHitRec[realRoots, 4];
SortHitRec[realRoots];
class.count ← 4;
class.params[1] ← realRoots[0];
class.params[2] ← realRoots[1];
class.params[3] ← realRoots[2];
class.params[4] ← realRoots[3];
class.surfaces[1] ← class.surfaces[2] ← class.surfaces[3] ← class.surfaces[4] ← surface;
class.classifs[1] ← FALSE; class.classifs[2] ← TRUE; class.classifs[3] ← FALSE;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.primitives[1] ← class.primitives[2] ← class.primitives[3] ← class.primitives[4] ← prim;
Calculate the surface normals at the four hits points.
x ← n + class.params[1]*N; y ← q + class.params[1]*Q; z ← l + class.params[1]*L;
theta ← RealFns.ArcTan[-z,x];
class.normals[1] ← [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z];
To normalize, divide by C. (I won't do this since the lighting model normalizes anyway)
x ← n + class.params[2]*N; y ← q + class.params[2]*Q; z ← l + class.params[2]*L;
theta ← RealFns.ArcTan[-z,x];
class.normals[2] ← [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z];
x ← n + class.params[3]*N; y ← q + class.params[3]*Q; z ← l + class.params[3]*L;
theta ← RealFns.ArcTan[-z,x];
class.normals[3] ← [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z];
x ← n + class.params[4]*N; y ← q + class.params[4]*Q; z ← l + class.params[4]*L;
theta ← RealFns.ArcTan[-z,x];
class.normals[4] ← [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z];
};
ENDCASE => ERROR;
};
SortHitRec: PROC [roots: REF Polynomial.RealRootRec] = {
temp: REAL;
FOR i: NAT IN[0..roots.nRoots-1) DO
FOR j: NAT IN[0..roots.nRoots-i-1) DO
IF roots[j] > roots[j+1] THEN {
temp ← roots[j];
roots[j] ← roots[j+1];
roots[j+1] ← temp;};
ENDLOOP;
ENDLOOP;
};
SortHitRec: PROC [roots: ARRAY[0..3] OF REAL, rootCount: NAT] = {
temp: REAL;
FOR i: NAT IN[0..rootCount-1) DO
FOR j: NAT IN[0..rootCount-i-1) DO
IF roots[j] > roots[j+1] THEN {
temp ← roots[j];
roots[j] ← roots[j+1];
roots[j+1] ← temp;};
ENDLOOP;
ENDLOOP;
};
Init: PROC = {
globalPoly ← Polynomial.Quartic[[0,0,0,0,0]];
};
Init[];
END.