[Cyan]<CedarChest6.0>SolidModeler>SteinerClassImpl.mesa!1

File: SteinerClassImpl.mesa

Last edited by Bier on January 10, 1985 5:22:13 pm PST

Author: Dennis Arnon and Eric Bier on January 10, 1985 10:07:59 pm PST

DIRECTORY

CastRays,

CSG,

DisplayList3d,

IO,

Polynomial,

Rope,

SV2d,

SV3d,

SVModelTypes,

SVRayTypes,

SVSceneTypes,

TFI3d;

SteinerClassImpl: PROGRAM

IMPORTS CastRays, CSG, DisplayList3d, IO, Polynomial, TFI3d

EXPORTS =

BEGIN

Camera: TYPE = SVModelTypes.Camera;

Point2d: TYPE = SV2d.Point2d;

Point3d: TYPE = SV3d.Point3d;

Shape: TYPE = SVSceneTypes.Shape;

Vector: TYPE = SV3d.Vector;

RAY CASTING TYPES

Assembly: TYPE = SVSceneTypes.Assembly;

Classification: TYPE = SVRayTypes.Classification;

Composite: TYPE = SVRayTypes.Composite;

MasterObject: TYPE = SVSceneTypes.MasterObject;

MasterObjectClass: TYPE = SVSceneTypes.MasterObjectClass;

MasterObjectClassList: TYPE = SVSceneTypes.MasterObjectClassList; -- LIST OF MasterObjectClass

Primitive: TYPE = SVRayTypes.Primitive;

Ray: TYPE = SVRayTypes.Ray;

globalPoly: Polynomial.Ref;

steinerClass: MasterObjectClass;

SteinerCast: PUBLIC PROC [cameraPoint: Point2d, localRay: Ray, masterObject: REF ANY, prim: Primitive]

RETURNS [class: Classification] = {

The Steiner surface is:

2 2 2 2 2 2
F(x,y,z) = y z + x z - 2 x y z + x y = 0

steinerData: SteinerRec;

realRoots: Polynomial.ShortRealRootRec; -- optional

realRoots: ARRAY[0..3] OF REAL;

realRootsOneOrigin: ARRAY[1..4] OF REAL;

rootCount: NAT;

x0, x1, y0, y1, z0, z1, x1y1: REAL;

x0x0, x1x1, y0y0, y1y1, z0z0, z1z1: REAL;

x0x1, y0y1, x0y0, x0y1, x1y0: REAL;

a4,a3,a2,a1,a0: REAL;

p: Point3d;

d: Vector;

SteinerGrad: PROC [x, y, z: REAL] RETURNS [Vector] ~ {

gradx, grady, gradz: REAL;

gradx ← 2. * ( x * z * z - y * z + x * y * y );

grady ← 2. * ( y * z * z - x * z + x * x * y );

gradz ← 2. * ( y * y * z + x * x * z - x * y );

RETURN [[gradx, grady, gradz]];

};

class ← CastRays.GetClassFromPool[];

steinerData ← NARROW[steinerRec];

[p, d] ← CSG.GetLocalRay[localRay]; -- p is base point, d is directiona

x0 ← p[1]; x1 ← d[1];

y0 ← p[2]; y1 ← d[2];

z0 ← p[3]; z1 ← d[3];

x0x0 ← x0 * x0;

x1x1 ← x1 * x1;

y0y0 ← y0 * y0;

y1y1 ← y1 * y1;

z0z0 ← z0 * z0;

z1z1 ← z1 * z1;

x0x1 ← x0 * x1;

y0y1 ← y0 * y1;

x0y0 ← x0 * y0;

x0y1 ← x0 * y1;

x1y0 ← x1 * y0;

The result of substituting

F, x: x1 * t + x0, y: y1 * t + y0, z: z1 * t + z0;

in Macsyma:

2 2 2 2 2 4
(d20)/R/ ((y1 + x1 ) z1 + x1 y1 ) t
2 2 2 2
+ ((2 y0 y1 + 2 x0 x1) z1 + ((2 y1 + 2 x1 ) z0 - 2 x1 y1) z1 + 2 x0 x1 y1
2 3 2 2 2
+ 2 x1 y0 y1) t + ((y0 + x0 ) z1 + ((4 y0 y1 + 4 x0 x1) z0 - 2 x0 y1
2 2 2 2 2
- 2 x1 y0) z1 + (y1 + x1 ) z0 - 2 x1 y1 z0 + x0 y1 + 4 x0 x1 y0 y1
2 2 2 2 2 2
+ x1 y0 ) t + (((2 y0 + 2 x0 ) z0 - 2 x0 y0) z1 + (2 y0 y1 + 2 x0 x1) z0
2 2 2 2 2
+ (- 2 x0 y1 - 2 x1 y0) z0 + 2 x0 y0 y1 + 2 x0 x1 y0 ) t + (y0 + x0 ) z0
2 2
- 2 x0 y0 z0 + x0 y0

a4 ← (y1y1 + x1x1 ) * z1z1 + x1x1 * y1y1 ;

a3 ← (2.*y0y1+2.*x0x1) * z1z1 +

((2.*y1y1 + 2.*x1x1)*z0 - 2.*x1y1)*z1 + 2.*x0x1*y1y1+ 2.*x1x1*y0y1 ;

a2 ← (y0y0 + x0x0 )*z1z1 + ((4.*y0y1 + 4.*x0x1)*z0 - 2.*x0y1 - 2.*x1y0)*z1

+ (y1y1 + x1x1 )*z0z0 - 2.*x1y1*z0 + x0x0*y1y1 + 4.*x0x1*y0y1 + x1x1 * y0y0;

a1 ← ((2.*y0y0 + 2.*x0x0)*z0 - 2.*x0y0)*z1 + (2.*y0y1 + 2.*x0x1)*z0z0

+ (-2.*x0y1 - 2.*x1y0)*z0 + 2.*x0x0*y0y1 + 2.*x0x1*y0y0;

a0 ← (y0y0 + x0x0 )*z0z0 - 2.*x0y0*z0 + x0x0*y0y0 ;

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 Steiner. Sort roots by size.

x, y, z: 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] ← NIL; class.surfaces[2] ← NIL;

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 hit points.

x ← x0 + class.params[1]*x1; y ← y0 + class.params[1]*y1; z ← z0 + class.params[1]*z1;

class.normals[1] ← SteinerGrad[x, y, z];

To normalize, divide by norm. (I won't do this since the lighting model normalizes anyway)

x ← x0 + class.params[2]*x1; y ← y0 + class.params[2]*y1; z ← z0 + class.params[2]*z1;

class.normals[2] ← SteinerGrad[x, y, z];

};

3 => {-- cuts one cross section. Skims another. For now, count as a miss.

class.count ← 0;

class.classifs[1] ← FALSE;

};

4 => {-- cuts two cross sections. Sort the roots. Pair them in order.

x, y, z: 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] ← NIL;

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 hit points.

x ← x0 + class.params[1]*x1; y ← y0 + class.params[1]*y1; z ← z0 + class.params[1]*z1;

class.normals[1] ← SteinerGrad[x, y, z];

To normalize, divide by norm. (I won't do this since the lighting model normalizes anyway)

x ← x0 + class.params[2]*x1; y ← y0 + class.params[2]*y1; z ← z0 + class.params[2]*z1;

class.normals[2] ← SteinerGrad[x, y, z];

x ← x0 + class.params[3]*x1; y ← y0 + class.params[3]*y1; z ← z0 + class.params[3]*z1;

class.normals[3] ← SteinerGrad[x, y, z];

To normalize, divide by norm. (I won't do this since the lighting model normalizes anyway)

x ← x0 + class.params[4]*x1; y ← y0 + class.params[4]*y1; z ← z0 + class.params[4]*z1;

class.normals[4] ← SteinerGrad[x, y, z];

};

ENDCASE => ERROR;

};

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 ← Polynomial.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;

};

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;

};

MakeMasterObject: PROC [name: Rope.ROPE] RETURNS [mo: MasterObject] = {

mainBody: REF ANY ← NIL;

lineBody: REF ANY ← NIL;

shadeBody: REF ANY ← lineBody;

rayCastBody: REF ANY ← NIL;

mo ← DisplayList3d.CreateMasterObject[name, steinerClass, mainBody, lineBody, shadeBody, rayCastBody];

};

SteinerFileout: PUBLIC PROC [f: IO.STREAM, mo: MasterObject] = {

Steiners can be built from scratch.

f.PutChar[IO.TAB];

f.PutF["data: procedural\n"];

};

SteinerFilein: PUBLIC PROC [f: IO.STREAM, name: Rope.ROPE] RETURNS [mo: MasterObject] = {

TFI3d.ReadRope[f, "data: procedural"];

TFI3d.ReadBlank[f];

mo ← MakeMasterObject[name];

};

Init: PROC = {

steiner: MasterObject;

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

steinerClass ← DisplayList3d.RegisterMasterObjectClass[

"steiner",

SteinerFilein,

SteinerFileout,

DisplayList3d.NoOpFileoutPoly,

SteinerCast,

DisplayList3d.NoOpRayCastNoBBoxes,

DisplayList3d.NoOpRayCastBoundingSpheres,

DisplayList3d.NoOpBoundHedron,

DisplayList3d.NoOpPreprocess,

DisplayList3d.NoOpLineDraw,

DisplayList3d.NoOpNormalsDraw,

DisplayList3d.NoOpCountPlanarSurfaces,

DisplayList3d.NoOpGetPlanarSurfaces,

DisplayList3d.NoOpDrawPlanarSurface,

DisplayList3d.NoOpDrawSubBoxes,

DisplayList3d.NoOpDrawSubSpheres];

steiner ← MakeMasterObject["steiner"];

DisplayList3d.RegisterMasterObject[steiner];

};

Init[];

END.