[Indigo]<Solidviews>SolidModeler>SVObjectCastImplC.mesa!4

File: SVObjectCastImplC.mesa

Last edited by Bier on January 2, 1983 8:31 pm

Author: Bier on June 1, 1987 11:45:31 am PDT

DIRECTORY

SVMasterObject, SVCastRays, SVRay, SVMasterObjectTypes, SVObjectCast, Polynomial, RealFns, SV2d, SV3d, SVSceneTypes, Roots;

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;

Debugging

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]

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;

};

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;

};

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

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;

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;

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.