DIRECTORY BasicObject3d, CastRays, CSG, ObjectCast, Polynomial, RealFns, SV2d, Roots; ObjectCastImplC: PROGRAM IMPORTS CastRays, Polynomial, RealFns EXPORTS ObjectCast = BEGIN 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: 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 _ Polynomial.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] = { 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; 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]; 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[a4, a3, a2, a1, a0]; 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, theta: 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; theta _ RealFns.ArcTan[-z,x]; class.normals[1] _ [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z]; 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]; 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[1] _ TRUE; class.classifs[2] _ 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; theta _ RealFns.ArcTan[-z,x]; class.normals[1] _ [-B*RealFns.Cos[theta] + x, y, B*RealFns.Sin[theta] + z]; 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: 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; }; Init: PROC = { globalPoly _ Polynomial.Quartic[[0,0,0,0,0]]; }; Init[]; END. lFile: ObjectCastImplC.mesa Last edited by Bier on January 2, 1983 8:31 pm Author: Bier on August 19, 1983 1:26 pm 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) 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; }; Êz– "Cedar" style˜Ihead1šœ™Jšœ.™.Jšœ(™(šÏk ˜ Jšœ˜J˜ Jšœ˜J˜Jšœ˜J˜Jšœ˜Jšœ˜—šœ˜Jšœ˜%Jšœ˜—Jšœ˜Jšœ€Ïb4œžeœÁž$œÃ™ÃJšœœœ˜*Jšœ œ˜Jšœœœ˜ Jšœœœ˜Jšœ œœ ˜Jšœ œ˜)Jšœ˜J˜šÏn œ œœœ-˜hJ™Iprocšœœ˜Lšœ ˜ Jšœ ˜ Jšœ ˜ Jšœ ˜ Jšœ ˜ Lšœ?˜?LšœH™HJšœK™KJšœL™LLšœœœ˜$šœ˜šœœ˜&šœœœ˜(Lšœ0˜0—Lšœ˜Lšœ(˜(L˜—Lšœ ˜—Lšœ˜Jšœ˜—š Ÿœœœ-œœ˜_Jšœ˜#J˜Jšœ4˜4J™J™(J™Jšœœœ˜ Jš œœœœœ˜Jšœœœ˜ Jšœœœœ˜Jšœœœœœœœœ˜;Jšœœ˜Jšœ$˜$Jšœœ˜Jšœ˜Jšœ˜Jšœœ˜2Jšœœ˜2Jšœœ˜2Jšœœœ˜ Jšœœœ˜ Jšœœœ˜Jšœœœ˜Jšœœœ˜ Jšœœœœœœœœœœ˜@Jšœœœœœœœœœ˜+Jšœœ œ œ˜Jšœœœœœœœœœ˜3Jšœœœœœœœœœ˜MJšœœœœœœœœœœœœœœœ˜JšœAœœœ˜bJšœœœœœœœœœ˜`J˜Jšœ.˜.J™Jšœ˜J™šœÏc˜Jšœ˜Jšœœ˜Jšœ˜—šœ ˜"Jšœ˜Jšœœ˜Jšœ˜—šœ ?˜EJ˜J˜šœ/˜5Jšœ)˜)Jšœ)˜)—š˜Jšœ)˜)Jšœ*˜*—Jšœ9˜9Jšœœœœ˜QJšœ1˜1J™5JšœR˜RJ˜JšœL˜LJ™XJšœR˜RJ˜JšœL˜LJšœ˜—šœ G˜MJšœ˜Jšœœ˜Jšœ˜—šœ B˜HJšœœ˜Jšœ™Jšœ˜J˜Jšœ(˜(Jšœ(˜(Jšœ(˜(Jšœ(˜(JšœX˜XJšœœœœ˜QJšœœœ˜5Jšœ]˜]J™6JšœR˜RJ˜JšœL˜LJ™XJšœR˜RJ˜JšœL˜LJšœR˜RJ˜JšœL˜LJšœR˜RJ˜JšœL˜LJ˜—Jšœœ˜Jšœ˜—šŸ œœ)˜9Jšœœ˜šœœœ˜#šœœœ˜%šœ)œ˜1Jšœ˜Jšœ(˜(Jšœ˜——Jšœ˜—Jšœ˜Jšœ˜—šŸ œœ1™AJšœœ™šœœœ™ šœœœ™"šœœ™Jšœ™Jšœ™Jšœ™——Jšœ™—Jšœ™Jšœ™—šŸœœ˜Jšœ-˜-Jšœ˜—J˜Jšœ˜Jšœ˜—…—ô+Ú