File: SVObjectCastImplC.mesa
Last edited by Bier on January 2, 1983 8:31 pm
Copyright © 1984 by Xerox Corporation. All rights reserved.
Author: Bier on August 29, 1987 2:10:52 pm PDT
DIRECTORY
SVMasterObject, SVCastRays, SVRay, SVMasterObjectTypes, SVObjectCast, Polynomial, RealFns, SV2d, SV3d, SVSceneTypes;
SVObjectCastImplC:
CEDAR PROGRAM
IMPORTS SVCastRays, SVRay, Polynomial, RealFns
EXPORTS SVObjectCast =
BEGIN
Classification: TYPE = SVSceneTypes.Classification;
Point2d: TYPE = SV2d.Point2d;
Point3d: TYPE = SV3d.Point3d;
Primitive: TYPE = SVSceneTypes.Primitive;
Ray: TYPE = SVSceneTypes.Ray;
Surface: TYPE = SVSceneTypes.Surface;
TorusRec: TYPE = SVMasterObjectTypes.TorusRec;
Vector3d: TYPE = SV3d.Vector3d;
ShortRealRootRec: TYPE = Polynomial.ShortRealRootRec;
Ref: TYPE = Polynomial.Ref;
globalPoly: Polynomial.Ref;
newWay: BOOL ← TRUE;
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:
PositiveRoots:
PROCEDURE [a4, a3, a2, a1, a0:
REAL]
RETURNS [rootArray: Polynomial.ShortRealRootRec] = {
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 ← CheapRealRoots[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.realRoot[i-1] <= 0
THEN {
FOR j:
NAT
IN [i..rootArray.nRoots-1]
DO
rootArray.realRoot[j-1] ← rootArray.realRoot[j];
ENDLOOP;
rootArray.nRoots ← rootArray.nRoots - 1;
}
ELSE {i ← i + 1};
ENDLOOP;
};
ToroidCast:
PUBLIC
PROC [localRay: Ray, prim: Primitive, torusRec:
REF
ANY, surface: Surface, positiveTOnly:
BOOL ←
TRUE]
RETURNS [class: Classification] = {
IF newWay THEN RETURN[NewToroidCast[localRay, prim, torusRec, surface]]
ELSE RETURN[OldToroidCast[localRay, prim, torusRec, surface]];
};
OldToroidCast:
PUBLIC
PROC [localRay: Ray, prim: Primitive, torusRec:
REF
ANY, surface: Surface]
RETURNS [class: Classification] = {
torusData: TorusRec;
realRoots: Polynomial.ShortRealRootRec; -- 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;
p: Point3d;
d: Vector3d;
class ← SVCastRays.GetClassFromPool[];
torusData ← NARROW[torusRec];
[p, d] ← SVRay.GetLocalRay[localRay];
B ← torusData.bigRadius;
C ← torusData.sectionRadius;
n ← p[1]; N ← d[1];
q ← p[2]; Q ← d[2];
l ← p[3]; L ← d[3];
BB ← B*B;
CC ← C*C;
D ← BB + CC;
E ← BB - CC;
F ← BB*CC;
nn ← n*n; NN ← N*N; qq ← q*q; QQ ← Q*Q; ll ← l*l; LL ← L*L;
NNNN ← NN*NN; QQQQ ← QQ*QQ; LLLL ← LL*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[1.0, a3/a4, a2/a4, a1/a4, a0/a4];
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, r, OneMinusBoverR: REAL;
class.count ← 2;
IF realRoots.realRoot[0] < realRoots.realRoot[1]
THEN
{class.params[1] ← realRoots.realRoot[0];
class.params[2] ← realRoots.realRoot[1]}
ELSE
{class.params[2] ← realRoots.realRoot[0];
class.params[1] ← realRoots.realRoot[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;
r ← RealFns.SqRt[x*x + z*z];
OneMinusBoverR ← 1 - (B/r);
class.normals[1] ← [x*OneMinusBoverR, y, z*OneMinusBoverR];
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;
r ← RealFns.SqRt[x*x + z*z];
OneMinusBoverR ← 1 - (B/r);
class.normals[2] ← [x*OneMinusBoverR, y, z*OneMinusBoverR];
};
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, r, OneMinusBoverR: REAL;
SortHitRec[realRoots, 4];
SortHitRec[realRoots];
class.count ← 4;
class.params[1] ← realRoots.realRoot[0];
class.params[2] ← realRoots.realRoot[1];
class.params[3] ← realRoots.realRoot[2];
class.params[4] ← realRoots.realRoot[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[4] ← TRUE; class.classifs[5] ← 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;
r ← RealFns.SqRt[x*x + z*z];
OneMinusBoverR ← 1 - (B/r);
class.normals[1] ← [x*OneMinusBoverR, y, z*OneMinusBoverR];
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;
r ← RealFns.SqRt[x*x + z*z];
OneMinusBoverR ← 1 - (B/r);
class.normals[2] ← [x*OneMinusBoverR, y, z*OneMinusBoverR];
x ← n + class.params[3]*N; y ← q + class.params[3]*Q; z ← l + class.params[3]*L;
r ← RealFns.SqRt[x*x + z*z];
OneMinusBoverR ← 1 - (B/r);
class.normals[3] ← [x*OneMinusBoverR, y, z*OneMinusBoverR];
x ← n + class.params[4]*N; y ← q + class.params[4]*Q; z ← l + class.params[4]*L;
r ← RealFns.SqRt[x*x + z*z];
OneMinusBoverR ← 1 - (B/r);
class.normals[4] ← [x*OneMinusBoverR, y, z*OneMinusBoverR];
};
ENDCASE => ERROR;
};
SortHitRec:
PROC [roots: Polynomial.ShortRealRootRec] = {
A simple bubble sort. Sorts in increasing order (depth).
temp: REAL;
FOR i:
NAT
IN[0..roots.nRoots-1)
DO
FOR j:
NAT
IN[0..roots.nRoots-i-1)
DO
IF roots.realRoot[j] > roots.realRoot[j+1]
THEN {
temp ← roots.realRoot[j];
roots.realRoot[j] ← roots.realRoot[j+1];
roots.realRoot[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;
};
EvalTorus:
PRIVATE
PROC [
RR, RRminusSS:
REAL, degree:
NAT, t:
REAL, localRay: Ray]
RETURNS [f:
REAL] = {
HH, x, y, z: REAL;
p1,p2,p3,d1,d2,d3: REAL;
basePt: Point3d;
direction: Vector3d;
[basePt, direction] ← SVRay.GetLocalRay[localRay];
p1 ← basePt[1];
p2 ← basePt[2];
p3 ← basePt[3];
d1 ← direction[1];
d2 ← direction[2];
d3 ← direction[3];
x ← p1 + t*d1;
y ← p2 + t*d2;
z ← p3 + t*d3;
SELECT degree FROM
4=>{
xxPLUSzz: REAL;
xxPLUSzz ← x*x + z*z;
HH ← xxPLUSzz + y*y;
f ← RRminusSS*RRminusSS + HH*(HH+2.0*RRminusSS) - 4.0*RR*xxPLUSzz;
};
3=>{
TwoHHdot: REAL;
HH ← x*x + y*y + z*z;
TwoHHdot ← 2.0*(x*d1 + y*d2 + z*d3);
f ← TwoHHdot*(HH + 2.0*RRminusSS) + HH*TwoHHdot - 8.0*RR*(x*d1 + z*d3);
};
2=>{
HHdot, HHdotPrime: REAL;
HH ← x*x + y*y + z*z;
HHdot ← (x*d1 + y*d2 + z*d3);
HHdotPrime ← (d1*d1+d2*d2+d3*d3);
f ← 4.0*HHdotPrime*(HH + RRminusSS) + 8.0*(HHdot*HHdot - RR*(d1*d1+d3*d3));
};
1=>{
f ← 24.0*(x*d1 + y*d2 + z*d3)*(d1*d1+d2*d2+d3*d3);
};
ENDCASE => ERROR;
};
TorusRoots:
PRIVATE
PROC [
RR, RRminusSS:
REAL, localRay: Ray, degree:
NAT, bound:
REAL]
RETURNS [roots: Polynomial.ShortRealRootRec] =
BEGIN
roots.nRoots ← 0;
IF degree=1
THEN
BEGIN
p: Point3d;
d: Vector3d;
[p,d] ← SVRay.GetLocalRay[localRay];
roots.nRoots ← 1;
roots.realRoot[0] ← -(p[1]*d[1] + p[2]*d[2] + p[3]*d[3])/(d[1]*d[1] + d[2]*d[2] + d[3]*d[3])
END
ELSE
BEGIN
BEGIN
extrema: Polynomial.ShortRealRootRec ← TorusRoots[RR, RRminusSS, localRay, degree-1, bound];
x: REAL;
debug Feedback.Append[IO.PutFR["%d degree, %d roots: ", IO.int[degree - 1], IO.int[extrema.nRoots]], TRUE, FALSE];
debug FOR i: NAT IN [0..extrema.nRoots) DO
debug Feedback.Append[IO.PutFR["%g ", IO.real[extrema.realRoot[i]]], TRUE, FALSE];
debug ENDLOOP;
debug Feedback.AppendTypescript[""];
IF extrema.nRoots>0
THEN
BEGIN
x ← TRootBetween[RR, RRminusSS, degree, -bound, extrema.realRoot[0], localRay];
IF x<=extrema.realRoot[0] THEN {roots.nRoots ← 1; roots.realRoot[0] ← x};
FOR i:
NAT
IN [0..extrema.nRoots-1)
DO
x ← TRootBetween[RR, RRminusSS, degree, extrema.realRoot[i], extrema.realRoot[i+1], localRay];
IF x<=extrema.realRoot[i+1] THEN {roots.realRoot[roots.nRoots] ← x; roots.nRoots ← roots.nRoots + 1};
ENDLOOP;
x ← TRootBetween[RR, RRminusSS, degree, extrema.realRoot[extrema.nRoots-1],bound, localRay];
IF x<=bound THEN {roots.realRoot[roots.nRoots] ← x; roots.nRoots ← roots.nRoots + 1}
END
ELSE
BEGIN
x ← TRootBetween[RR, RRminusSS, degree, -bound, bound, localRay];
IF x<=bound THEN {roots.realRoot[0] ← x; roots.nRoots ← 1}
END
END
END
END;
TRootBound:
PROC [
R,
S:
REAL, p: Point3d, d: Vector3d]
RETURNS [
REAL] =
BEGIN
-- finds a bound for the absolute values of the roots
We do this for the torus by reasoning geometrically. In local coordinates, the torus is
centered on [0,0,0]. The ray base-point is at point p. Hence (|p| + R + S)/|d| will bound
the value of d. We can use the 1-norm to overestimate |p| and the sup-norm to underestimate |d| to get a conservative bound on t.
p1, dSup: REAL;
p1 ← ABS[p[1]]+ABS[p[2]]+ABS[p[3]];
dSup ← MAX[ABS[d[1]], MAX[ABS[d[2]], ABS[d[3]]]];
RETURN[(p1 + R + S)/dSup];
END;
TRootBetween:
PROC [
RR, RRminusSS:
REAL, degree:
NAT, x0, x1:
REAL, localRay: Ray]
RETURNS [x: REAL] =
BEGIN
y0: REAL ← EvalTorus[RR, RRminusSS, degree, x0, localRay];
y1: REAL ← EvalTorus[RR, RRminusSS, degree, x1, localRay];
debug xi: ARRAY [0..50) OF REAL;
debug i: NAT ← 0;
xx: REAL;
xx0: REAL ← IF y0<y1 THEN x0 ELSE x1;
xx1: REAL ← IF y0<y1 THEN x1 ELSE x0;
yy0: REAL ← MIN[y0,y1];
yy1: REAL ← MAX[y0,y1];
NewtonStep:
PROC [x:
REAL]
RETURNS [newx:
REAL] =
BEGIN
y: REAL ← EvalTorus[RR, RRminusSS, degree, x, localRay];
deriv: REAL ← EvalTorus[RR, RRminusSS, degree-1, x, localRay];
IF y<0 THEN {IF y>yy0 THEN {xx0 ← x; yy0 ← y}};
IF y>0 THEN {IF y<yy1 THEN {xx1 ← x; yy1 ← y}};
newx ← IF deriv=0 THEN x+y ELSE x-y/deriv;
END;
IF y0=0 THEN RETURN[x0];
IF y1=0 THEN RETURN[x1];
IF (y0 > 0.0 AND y1 > 0.0) OR (y0 < 0.0 AND y1 < 0.0) THEN RETURN[x1+100]; -- return a non-sense value which the caller will check for (there is no root in this interval).
IF y0*y1 > 0 THEN RETURN[x1+100]; -- return a non-sense value which the caller will check for (there is no root in this interval).
xx ← x ← (x0+x1)/2;
first try Newton's method
WHILE x
IN [x0..x1]
DO
newx:REAL ← NewtonStep[x];
IF x=newx
THEN {
debug Feedback.Append[IO.PutFR["TNewton: %d steps: ", IO.int[i]], TRUE,FALSE];
debug FOR j: NAT IN [0..i) DO
debug Feedback.Append[IO.PutFR["%g ", IO.real[xi[j]]], TRUE, FALSE];
debug ENDLOOP;
debug Feedback.AppendTypescript[""];
RETURN;
};
x ← NewtonStep[newx];
xx ← NewtonStep[xx];
debug xi[i] ← xx; i ← i+1;
IF xx=x THEN EXIT;
ENDLOOP;
If it oscillated or went out of range, try bisection
debug SIGNAL Bisection;
IF xx0
IN [x0..x1]
AND xx1
IN [x0..x1]
THEN
BEGIN
x0 ← MIN[xx0,xx1];
x1 ← MAX[xx0,xx1];
END;
y0 ← EvalTorus[RR, RRminusSS, degree, x0, localRay];
y1 ← EvalTorus[RR, RRminusSS, degree, x1, localRay];
THROUGH [0..500)
DO
y: REAL ← EvalTorus[RR, RRminusSS, degree, x ← (x0+x1)/2, localRay];
IF x=x0
OR x=x1
THEN
{IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]};
IF y*y0 < 0 THEN {x1 ← x; y1 ← y}
ELSE {x0 ← x; y0 ← y};
IF y0*y1>=0 THEN {IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]}
ENDLOOP;
ERROR;
END;
PositiveTorusRoots:
PRIVATE
PROC [
RR, RRminusSS:
REAL, localRay: Ray, bound:
REAL]
RETURNS [roots: Polynomial.ShortRealRootRec] = {
i: NAT ← 1;
roots ← TorusRoots[RR, RRminusSS, localRay, 4, bound];
debug Feedback.Append[IO.PutFR["%d degree, %d roots: ", IO.int[4], IO.int[roots.nRoots]], TRUE, FALSE];
debug FOR i: NAT IN [0..roots.nRoots) DO
debug Feedback.Append[IO.PutFR["%g ", IO.real[roots.realRoot[i]]], TRUE, FALSE];
debug ENDLOOP;
debug Feedback.AppendTypescript[""];
IF roots.nRoots = 0 THEN RETURN;
WHILE i <= roots.nRoots
DO
IF roots.realRoot[i-1] <= 0
THEN {
FOR j:
NAT
IN [i..roots.nRoots-1]
DO
roots.realRoot[j-1] ← roots.realRoot[j];
ENDLOOP;
roots.nRoots ← roots.nRoots - 1;
}
ELSE {i ← i + 1};
ENDLOOP;
};
InsideTorus:
PRIVATE
PROC [
R,
S:
REAL, testPt: Point3d]
RETURNS [
BOOL] = {
r, rMinusR: REAL;
r ← RealFns.SqRt[testPt[1]*testPt[1] + testPt[3]*testPt[3]];
rMinusR ← r-R;
RETURN[rMinusR*rMinusR + testPt[2]*testPt[2] <= S*S];
};
NewToroidCast:
PUBLIC
PROC [localRay: Ray, prim: Primitive, torusRec:
REF
ANY, surface: Surface]
RETURNS [class: Classification] = {
torusData: TorusRec;
realRoots: Polynomial.ShortRealRootRec; -- optional
B, R, S, RR, RRminusSS: REAL; -- big and little torus radii
p: Point3d;
d: Vector3d;
bound: REAL;
class ← SVCastRays.GetClassFromPool[];
torusData ← NARROW[torusRec];
R ← B ← torusData.bigRadius;
S ← torusData.sectionRadius;
RR ← R*R;
RRminusSS ← RR - S*S;
[p, d] ← SVRay.GetLocalRay[localRay];
bound ← TRootBound[R, S, p, d];
realRoots ← PositiveTorusRoots[RR, RRminusSS, localRay, bound];
SELECT realRoots.nRoots FROM
SELECT rootCount FROM
0 => {
-- complete miss.
class.count ← 0;
class.classifs[1] ← FALSE;
};
1 => {
-- a skim. Count as a miss.
IF InsideTorus[
R,
S, p]
THEN {
x, y, z, r, OneMinusBoverR: REAL;
class.count ← 1;
class.classifs[1] ← TRUE; class.classifs[2] ← FALSE;
class.params[1] ← realRoots.realRoot[0];
Calculate the surface normals at the two hits points.
[x, y, z] ← SVRay.EvaluateLocalRay2[localRay, class.params[1]];
r ← RealFns.SqRt[x*x + z*z];
OneMinusBoverR ← 1 - (B/r);
class.normals[1] ← [x*OneMinusBoverR, y, z*OneMinusBoverR];
class.surfaces[1] ← surface;
class.primitives[1] ← prim;
}
ELSE SVCastRays.MakeClassAMiss[class];
};
2 => {
-- cuts one cross section of the doughnut. Sort roots by size.
x, y, z, r, OneMinusBoverR: REAL;
class.count ← 2;
IF realRoots.realRoot[0] < realRoots.realRoot[1]
THEN
{class.params[1] ← realRoots.realRoot[0];
class.params[2] ← realRoots.realRoot[1]}
ELSE
{class.params[2] ← realRoots.realRoot[0];
class.params[1] ← realRoots.realRoot[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, y, z] ← SVRay.EvaluateLocalRay2[localRay, class.params[1]];
r ← RealFns.SqRt[x*x + z*z];
OneMinusBoverR ← 1 - (B/r);
class.normals[1] ← [x*OneMinusBoverR, y, z*OneMinusBoverR];
To normalize, divide by C. (I won't do this since the lighting model normalizes anyway)
[x, y, z] ← SVRay.EvaluateLocalRay2[localRay, class.params[2]];
r ← RealFns.SqRt[x*x + z*z];
OneMinusBoverR ← 1 - (B/r);
class.normals[2] ← [x*OneMinusBoverR, y, z*OneMinusBoverR];
};
3 => {
-- cuts one cross section. Skims the other (or starts inside of the torus). 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, r, OneMinusBoverR: REAL;
SortHitRec[realRoots, 4];
SortHitRec[realRoots];
class.count ← 4;
class.params[1] ← realRoots.realRoot[0];
class.params[2] ← realRoots.realRoot[1];
class.params[3] ← realRoots.realRoot[2];
class.params[4] ← realRoots.realRoot[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[4] ← TRUE; class.classifs[5] ← FALSE;
class.primitives[1] ← class.primitives[2] ← class.primitives[3] ← class.primitives[4] ← prim;
Calculate the surface normals at the four hits points.
[x, y, z] ← SVRay.EvaluateLocalRay2[localRay, class.params[1]];
r ← RealFns.SqRt[x*x + z*z];
OneMinusBoverR ← 1 - (B/r);
class.normals[1] ← [x*OneMinusBoverR, y, z*OneMinusBoverR];
To normalize, divide by C. (I won't do this since the lighting model normalizes anyway)
[x, y, z] ← SVRay.EvaluateLocalRay2[localRay, class.params[2]];
r ← RealFns.SqRt[x*x + z*z];
OneMinusBoverR ← 1 - (B/r);
class.normals[2] ← [x*OneMinusBoverR, y, z*OneMinusBoverR];
[x, y, z] ← SVRay.EvaluateLocalRay2[localRay, class.params[3]];
r ← RealFns.SqRt[x*x + z*z];
OneMinusBoverR ← 1 - (B/r);
class.normals[3] ← [x*OneMinusBoverR, y, z*OneMinusBoverR];
[x, y, z] ← SVRay.EvaluateLocalRay2[localRay, class.params[4]];
r ← RealFns.SqRt[x*x + z*z];
OneMinusBoverR ← 1 - (B/r);
class.normals[4] ← [x*OneMinusBoverR, y, z*OneMinusBoverR];
};
ENDCASE => ERROR;
};
RootBound:
PROC [poly: Ref]
RETURNS [
REAL] =
BEGIN -- finds a bound for the absolute values of the roots
d: NAT ← Degree[poly];
b: REAL ← 0;
FOR i:NAT IN [0..d) DO b ← b + ABS[poly[i]] ENDLOOP;
RETURN [b*1.01/ABS[poly[d]]+1];
END;
RootBetween:
PROC [poly: Ref, x0, x1:
REAL]
RETURNS [x: REAL] =
BEGIN
debug xi: ARRAY [0..50) OF REAL;
debug i: NAT ← 0;
y0: REAL ← Eval[poly,x0];
y1: REAL ← Eval[poly,x1];
xx: REAL;
xx0: REAL ← IF y0<y1 THEN x0 ELSE x1;
xx1: REAL ← IF y0<y1 THEN x1 ELSE x0;
yy0: REAL ← MIN[y0,y1];
yy1: REAL ← MAX[y0,y1];
NewtonStep:
PROC [x:
REAL]
RETURNS [newx:
REAL] =
BEGIN
y: REAL ← Eval[poly, x];
deriv: REAL ← EvalDeriv[poly, x];
IF y<0 THEN {IF y>yy0 THEN {xx0 ← x; yy0 ← y}};
IF y>0 THEN {IF y<yy1 THEN {xx1 ← x; yy1 ← y}};
newx ← IF deriv=0 THEN x+y ELSE x-y/deriv;
END;
IF y0=0 THEN RETURN[x0];
IF y1=0 THEN RETURN[x1];
IF (y0 > 0.0 AND y1 > 0.0) OR (y0 < 0.0 AND y1 < 0.0) THEN RETURN[x1+100]; -- return a non-sense value which the caller will check for (there is no root in this interval).
IF y0*y1 > 0 THEN RETURN[x1+100];
xx ← x ← (x0+x1)/2;
first try Newton's method
WHILE x
IN [x0..x1]
DO
newx:REAL ← NewtonStep[x];
IF x=newx
THEN
{
debug Feedback.Append[IO.PutFR["Newton: %d steps: ", IO.int[i]], TRUE,FALSE];
debug FOR j: NAT IN [0..i) DO
debug Feedback.Append[IO.PutFR["%g ", IO.real[xi[j]]], TRUE, FALSE];
debug ENDLOOP;
debug Feedback.AppendTypescript[""];
RETURN;
};
x ← NewtonStep[newx];
xx ← NewtonStep[xx];
debug xi[i] ← xx; i ← i+1;
IF xx=x THEN EXIT;
ENDLOOP;
If it oscillated or went out of range, try bisection
debug SIGNAL Bisection;
IF xx0
IN [x0..x1]
AND xx1
IN [x0..x1]
THEN
BEGIN
x0 ← MIN[xx0,xx1];
x1 ← MAX[xx0,xx1];
END;
y0 ← Eval[poly,x0];
y1 ← Eval[poly,x1];
THROUGH [0..500)
DO
y: REAL ← Eval[poly, x ← (x0+x1)/2];
IF x=x0
OR x=x1
THEN
{IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]};
IF y*y0 < 0 THEN {x1 ← x; y1 ← y}
ELSE {x0 ← x; y0 ← y};
IF y0*y1>=0 THEN {IF ABS[y0] < ABS[y1] THEN RETURN[x0] ELSE RETURN[x1]}
ENDLOOP;
ERROR;
END;
Degree:
PUBLIC
PROC [poly: Ref]
RETURNS [
NAT] =
BEGIN
FOR i:
NAT
DECREASING
IN [0..poly.maxDegreePlusOne)
DO
IF poly[i] # 0 THEN RETURN[i];
ENDLOOP;
RETURN[0];
END;
Differentiate:
PUBLIC
PROC [poly: Ref] =
BEGIN
d:NAT ← Degree[poly];
FOR i: NAT IN (0..d] DO poly[i-1] ← i*poly[i] ENDLOOP;
poly[d] ← 0.0;
END;
CheapRealRoots:
PUBLIC
PROC [poly: Ref]
RETURNS [roots: ShortRealRootRec] =
BEGIN
d: NAT ← Degree[poly];
roots.nRoots ← 0;
IF d<=1
THEN
BEGIN
IF d=1 THEN {roots.nRoots ← 1; roots.realRoot[0] ← -poly[0]/poly[1]}
END
ELSE
BEGIN
bound: REAL ← RootBound[poly];
savedCoeff: ARRAY [0..5] OF REAL;
FOR i: NAT IN [0..d] DO savedCoeff[i] ← poly[i] ENDLOOP;
Differentiate[poly];
BEGIN
extrema: ShortRealRootRec ← CheapRealRoots[poly];
x: REAL;
debug Feedback.Append[IO.PutFR["%g bound, %d degree %d roots: ", IO.real[bound], IO.int[d-1], IO.int[extrema.nRoots]], TRUE, FALSE];
debug FOR i: NAT IN [0..extrema.nRoots) DO
debug Feedback.Append[IO.PutFR["%g ", IO.real[extrema.realRoot[i]]], TRUE, FALSE];
debug ENDLOOP;
debug Feedback.AppendTypescript[""];
FOR i: NAT IN [0..d] DO poly[i] ← savedCoeff[i] ENDLOOP;
IF extrema.nRoots>0
THEN
BEGIN
x ← RootBetween[poly,-bound,extrema.realRoot[0]];
IF x<=extrema.realRoot[0] THEN {roots.nRoots ← 1; roots.realRoot[0] ← x};
FOR i:
NAT
IN [0..extrema.nRoots-1)
DO
x ← RootBetween[poly,extrema.realRoot[i],extrema.realRoot[i+1]];
IF x<=extrema.realRoot[i+1] THEN {roots.realRoot[roots.nRoots] ← x; roots.nRoots ← roots.nRoots + 1};
ENDLOOP;
x ← RootBetween[poly,extrema.realRoot[extrema.nRoots-1],bound];
IF x<=bound THEN {roots.realRoot[roots.nRoots] ← x; roots.nRoots ← roots.nRoots + 1}
END
ELSE
BEGIN
x ← RootBetween[poly,-bound,bound];
IF x<=bound THEN {roots.realRoot[0] ← x; roots.nRoots ← 1}
END
END
END
END;
EvalDeriv:
PUBLIC
PROC [f: Ref, x:
REAL]
RETURNS [y:
REAL] =
BEGIN
y ← 0.0;
FOR i:
NAT
DECREASING
IN (0..Degree[f]]
DO
y ← i*f[i] + y*x
ENDLOOP
END;
Eval:
PUBLIC
PROC [f: Ref, x:
REAL]
RETURNS [y:
REAL] =
BEGIN
y ← 0.0;
FOR i:
NAT
DECREASING
IN [0..Degree[f]]
DO
y ← y*x + f[i]
ENDLOOP
END;
Init:
PROC = {
globalPoly ← Polynomial.Quartic[[0,0,0,0,0]];
};
Init[];
END.