DIRECTORY BasicObject3d, CastRays, CSG, ObjectCast, Polynomial, RealFns, SV2d, SV3d, SVRayTypes, Roots; ObjectCastImplC: PROGRAM IMPORTS CastRays, CSG, Polynomial, RealFns EXPORTS ObjectCast = BEGIN Classification: TYPE = SVRayTypes.Classification; Point2d: TYPE = SV2d.Point2d; Point3d: TYPE = SV3d.Point3d; Primitive: TYPE = SVRayTypes.Primitive; Ray: TYPE = SVRayTypes.Ray; Surface: TYPE = SVRayTypes.Surface; TorusRec: TYPE = BasicObject3d.TorusRec; Vector: TYPE = SV3d.Vector; ShortRealRootRec: TYPE = Polynomial.ShortRealRootRec; Ref: TYPE = Polynomial.Ref; globalPoly: Polynomial.Ref; newWay: BOOL _ TRUE; PositiveRoots: PROCEDURE [a4, a3, a2, a1, a0: REAL] RETURNS [rootArray: Polynomial.ShortRealRootRec] = { 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 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 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: Vector; class _ CastRays.GetClassFromPool[]; torusData _ NARROW[torusRec]; [p, d] _ CSG.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 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; 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]; 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]; 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; 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]; 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] = { 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; }; 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: Vector; [basePt, direction] _ CSG.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: Vector; [p,d] _ CSG.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; 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: Vector] RETURNS [REAL] = BEGIN -- finds a bound for the absolute values of the roots 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]; xx: REAL; xx0: REAL _ IF y0yy0 THEN {xx0 _ x; yy0 _ y}}; IF y>0 THEN {IF y 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). xx _ x _ (x0+x1)/2; WHILE x IN [x0..x1] DO newx:REAL _ NewtonStep[x]; IF x=newx THEN { RETURN; }; x _ NewtonStep[newx]; xx _ NewtonStep[xx]; IF xx=x THEN EXIT; ENDLOOP; 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]; 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: Vector; bound: REAL; class _ CastRays.GetClassFromPool[]; torusData _ NARROW[torusRec]; R _ B _ torusData.bigRadius; S _ torusData.sectionRadius; RR _ R*R; RRminusSS _ RR - S*S; [p, d] _ CSG.GetLocalRay[localRay]; bound _ TRootBound[R, S, p, d]; realRoots _ PositiveTorusRoots[RR, RRminusSS, localRay, bound]; SELECT realRoots.nRoots 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]; [x, y, z] _ CSG.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 CastRays.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; [x, y, z] _ CSG.EvaluateLocalRay2[localRay, class.params[1]]; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[1] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; [x, y, z] _ CSG.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]; 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; [x, y, z] _ CSG.EvaluateLocalRay2[localRay, class.params[1]]; r _ RealFns.SqRt[x*x + z*z]; OneMinusBoverR _ 1 - (B/r); class.normals[1] _ [x*OneMinusBoverR, y, z*OneMinusBoverR]; [x, y, z] _ CSG.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] _ CSG.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] _ CSG.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 y0: REAL _ Eval[poly,x0]; y1: REAL _ Eval[poly,x1]; xx: REAL; xx0: REAL _ IF y0yy0 THEN {xx0 _ x; yy0 _ y}}; IF y>0 THEN {IF y 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). xx _ x _ (x0+x1)/2; WHILE x IN [x0..x1] DO newx:REAL _ NewtonStep[x]; IF x=newx THEN { RETURN; }; x _ NewtonStep[newx]; xx _ NewtonStep[xx]; IF xx=x THEN EXIT; ENDLOOP; 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; 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. File: ObjectCastImplC.mesa Last edited by Bier on January 2, 1983 8:31 pm Copyright c 1984 by Xerox Corporation. All rights reserved. Author: Bier on January 2, 1985 1:13:56 pm PST Debugging 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: Use only the positive roots. [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]; realRoots: ARRAY[0..3] OF REAL; realRootsOneOrigin: ARRAY[1..4] OF REAL; rootCount: NAT; SELECT rootCount FROM Calculate the surface normals at the two hits points. To normalize, divide by C. (I won't do this since the lighting model normalizes anyway) SortHitRec[realRoots, 4]; Calculate the surface normals at the four hits points. To normalize, divide by C. (I won't do this since the lighting model normalizes anyway) A simple bubble sort. Sorts in increasing order (depth). 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; }; debug SVError.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 SVError.Append[IO.PutFR["%g ", IO.real[extrema.realRoot[i]]], TRUE, FALSE]; debug ENDLOOP; debug SVError.AppendTypescript[""]; 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. debug xi: ARRAY [0..50) OF REAL; debug i: NAT _ 0; 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). first try Newton's method debug SVError.Append[IO.PutFR["TNewton: %d steps: ", IO.int[i]], TRUE,FALSE]; debug FOR j: NAT IN [0..i) DO debug SVError.Append[IO.PutFR["%g ", IO.real[xi[j]]], TRUE, FALSE]; debug ENDLOOP; debug SVError.AppendTypescript[""]; debug xi[i] _ xx; i _ i+1; If it oscillated or went out of range, try bisection debug SIGNAL Bisection; debug SVError.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 SVError.Append[IO.PutFR["%g ", IO.real[roots.realRoot[i]]], TRUE, FALSE]; debug ENDLOOP; debug SVError.AppendTypescript[""]; SELECT rootCount FROM Calculate the surface normals at the two hits points. Calculate the surface normals at the two hits points. To normalize, divide by C. (I won't do this since the lighting model normalizes anyway) SortHitRec[realRoots, 4]; Calculate the surface normals at the four hits points. To normalize, divide by C. (I won't do this since the lighting model normalizes anyway) debug xi: ARRAY [0..50) OF REAL; debug i: NAT _ 0; IF y0*y1 > 0 THEN RETURN[x1+100]; first try Newton's method debug SVError.Append[IO.PutFR["Newton: %d steps: ", IO.int[i]], TRUE,FALSE]; debug FOR j: NAT IN [0..i) DO debug SVError.Append[IO.PutFR["%g ", IO.real[xi[j]]], TRUE, FALSE]; debug ENDLOOP; debug SVError.AppendTypescript[""]; debug xi[i] _ xx; i _ i+1; If it oscillated or went out of range, try bisection debug SIGNAL Bisection; debug SVError.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 SVError.Append[IO.PutFR["%g ", IO.real[extrema.realRoot[i]]], TRUE, FALSE]; debug ENDLOOP; debug SVError.AppendTypescript[""]; ΚP– "Cedar" style˜Ihead1šœ™Jšœ.™.Jšœ Οmœ1™J˜J˜——š   œžœžœ-žœžœ˜aJšžœ˜#J˜Jšœ)‘ ˜4J™J™(J™Jšžœžœžœ˜ Jš œžœžœžœžœ˜Jšžœžœžœ˜ Jšžœžœžœžœ˜ Jšœžœžœžœžœžœžœžœ˜;Jšœžœ˜J˜ J˜ Jšœ$˜$Jšœ žœ ˜Jšœ žœ˜#Jšžœ˜Jšžœ˜Jšœ žœ˜Jšœ žœ˜Jšœ žœ˜Jšžœžœžœ˜ Jšžœžœžœ˜ Jšžœžœžœ˜ Jšžœžœžœ˜ Jšžœžœžœ˜ Jšœ žœžœžœžœžœžœžœžœžœ˜@Jšžœžœžœžœžœžœžœžœžœ˜+Jšœžœ žœ žœ˜Jšœžœžœžœžœžœžœžœžœ˜3Jšœ žœžœžœžœžœžœžœžœ˜MJšœžœžœžœžœžœžœžœžœžœžœžœžœžœžœ˜JšœAžœžœžœ˜bJšœžœžœžœžœžœžœžœžœ˜`J˜Jšœ;˜;J™Jšžœž˜J™šœ‘˜Jšœ˜Jšœžœ˜Jšœ˜—šœ‘˜"Jšœ˜Jšœžœ˜Jšœ˜—šœ‘?˜EJšœžœ˜!J˜šžœ/ž˜5Jšœ)˜)Jšœ)˜)—šž˜Jšœ)˜)Jšœ*˜*—Jšœ9˜9Jšœžœžœžœ˜QJšœ1˜1J™5Jšœžœžœžœ˜RJ˜Jšœ˜Jšœ;˜;J™XJšœžœžœžœ˜RJ˜Jšœ˜Jšœ;˜;Jšœ˜—šœ‘G˜MJšœ˜Jšœžœ˜Jšœ˜—šœ‘B˜HJšœžœ˜!Jšœ™Jšœ˜J˜Jšœ(˜(Jšœ(˜(Jšœ(˜(Jšœ(˜(JšœX˜XJšœžœžœžœ˜QJšœžœžœ˜5Jšœ]˜]J™6Jšœžœžœžœ˜RJ˜Jšœ˜Jšœ;˜;J™XJšœžœžœžœ˜RJ˜Jšœ˜Jšœ;˜;Jšœžœžœžœ˜RJ˜Jšœ˜Jšœ;˜;Jšœžœžœžœ˜RJ˜Jšœ˜Jšœ;˜;J˜—Jšžœžœ˜Jšœ˜—J˜š  œžœ)˜9Jšœ8™8Jšœžœ˜ šžœžœžœž˜#šžœžœžœž˜%šžœ)žœ˜1Jšœ˜Jšœ(˜(Jšœ˜——Jšžœ˜—Jšžœ˜Jšœ˜—š  œžœ1™AJšœžœ™ šžœžœžœž™ šžœžœžœž™"šžœžœ™Jšœ™Jšœ™Jšœ™——Jšžœ™—Jšžœ™Jšœ™—š  œžœžœžœ žœ žœžœžœžœ˜hJšžœ žœ˜J˜J˜J˜Jšœ  œ ˜0Lšœ˜Lšœ˜Lšœ˜Lšœ˜Lšœ˜Lšœ˜Lšœ˜L˜L˜Jšžœž˜šœ˜Jšœ žœ˜J˜Jšžœ˜Jšœžœžœžœ ˜BJšœ˜—šœ˜Jšœ žœ˜Jšžœ˜J˜$Jšœžœžœžœ˜GJ˜—šœ˜Jšœžœ˜Jšžœ˜J˜Jšœ!˜!Jšœžœ#žœ˜KJšœ˜—šœ˜Jšœ2˜2J˜—Jšžœžœ˜J˜J˜—š  œžœžœžœ žœžœ žœžœ'˜†Jšž˜J˜šžœ ž˜šž˜J˜ Jšœ ˜ Jšœžœ˜"Jšœ˜Jšœ\˜\—Jšž˜—šž˜šž˜Jšž˜Jšœ2žœ(˜\Jšœžœ˜Jš œžœ žœžœžœžœ™qšœžœžœžœž™*Jš œžœžœžœžœ™Q—Jšœžœ™Jšœ#™#šžœž˜Jšž˜Jšœžœ<˜OJšžœžœ+˜Išžœžœžœž˜%JšœžœK˜^JšžœžœE˜eJšžœ˜—JšœžœI˜\Jšžœ žœD˜TJšž˜—šž˜Jšž˜Jšœžœ.˜AJšžœ žœ*˜:Jšž˜—Jšž˜—Jšž˜—Jšžœ˜J˜—J˜š  œžœžœžœžœžœžœ˜Ešžœ‘5˜;J™·Jšœ žœ˜J˜#Jšœžœ žœ˜1Jšžœžœžœ˜—Jšžœ˜J˜—š   œžœžœ žœ žœ žœ˜RJšžœžœ˜Jšž˜Jšœžœ žœ#˜:Jšœžœ žœ#˜:Jšœ!™!Jšœ™Jšœžœ˜ Jš œžœžœžœžœ˜%Jš œžœžœžœžœ˜%Jšœžœžœ˜Jšœžœžœ˜š   œžœžœžœžœ˜0Jšž˜Jšœžœ žœ"˜8Jšœžœ žœ$˜>Jšžœžœžœžœ˜/Jšžœžœžœžœ˜/Jšœžœ žœžœ ˜*Jšžœ˜—Jšžœžœžœ˜Jšžœžœžœ˜Jš žœ žœ žœ žœ žœžœ ‘`˜«Jšžœ žœžœ ‘`™‚J˜Jšœ™šžœžœ ž˜Jšœžœ˜šžœž˜Jš œžœžœ žœžœ™Nšœžœžœžœž™Jš œžœžœžœžœ™D—Jšœžœ™J™$Jšžœ˜J˜—J˜J˜Jšœ™Jšžœžœžœ˜Jšžœ˜—Jšœ4™4Jšœ™š žœžœ žœžœ ž˜+Jšž˜Jšœžœ ˜Jšœžœ ˜Jšžœ˜—Jšœžœ#˜4Jšœžœ#˜4šžœ ž˜Jšœžœ žœ.˜Dšžœžœž˜Jšœžœžœžœžœžœžœžœ˜7—Jšžœ žœ˜!Jšžœ˜Jšžœ žœžœžœžœžœžœžœžœ˜GJšžœ˜—Jšžœ˜Jšžœ˜J˜J˜—š œžœžœžœ žœžœžœ)˜ƒJšœžœ˜ Jšœžœ!˜6Jš œžœ žœ žœžœžœ™gšœžœžœžœž™)Jš œžœžœžœžœ™P—Jšœžœ™Jšœ$™$Jšžœžœžœ˜ šžœž˜šžœžœ˜"šžœžœžœž˜$Lšœ(˜(—Lšžœ˜Lšœ ˜ L˜—Lšžœ ˜—Lšžœ˜L˜L˜—š  œžœžœžœžœžœžœžœ˜JJšœ žœ˜J˜