[Cyan]<ISL>SolidViews>ObjectCastImplC.mesa!1

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];

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

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;

};

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;

};

Init: PROC = {

globalPoly ← Polynomial.Quartic[[0,0,0,0,0]];

};

Init[];

END.