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;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] = {realRoots: Polynomial.ShortRealRootRec;  -- optional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[];[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;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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;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];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];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;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];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];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] = {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;};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;};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] = {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.   „  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RAY CASTING TYPES 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: ARRAY[0..3] OF REAL;realRootsOneOrigin: ARRAY[1..4] OF REAL;rootCount: NAT;steinerData _ NARROW[steinerRec];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SELECT rootCount FROMCalculate the surface normals at the two hit points.To normalize, divide by norm.  (I won't do this since the lighting model normalizes anyway)SortHitRec[realRoots, 4];Calculate the surface normals at the four hit points.To normalize, divide by norm.  (I won't do this since the lighting model normalizes anyway)To normalize, divide by norm.  (I won't do this since the lighting model normalizes anyway)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];A simple bubble sort. Sorts in increasing order (depth).Steiners can be built from scratch. Ê4  ˜ Ihead1™J™6J™FJ˜ šÏk	˜	Jšœ˜	J˜Jšœ˜J˜J˜J˜J˜J˜Jšœ˜Jšœ˜Jšœ˜Jšœ˜J˜ —šœ˜Jšœœœ˜;Jšœ˜
—Jš˜J˜ Iprocšœœ˜#Lšœ	œ˜Lšœ	œ˜L˜ Lšœœ˜!Lšœœ˜L˜ Lšœ™˜ Lšœ
œ˜'Lšœœ˜2Lšœœ˜'Lšœœ˜/Lšœœ"˜9Lšœœ'Ïc˜^Lšœœ˜'Lšœœ˜L˜ Jšœ˜Jšœ ˜ J˜ —J˜ š
Ïnœœœ5œœ˜fšœ˜#J˜ JšÏf™Jšœ v™w—J™Jšœ)ž˜4J™J™(J™J˜#Jšœ$œ˜)Jš	œ œ œ œœ˜#Jšœœ˜J˜J˜
šŸœœœœ˜6Jšœœ˜J˜1J˜/J˜/Jšœ˜J˜—J˜ Jšœ$˜$Jšœœ™!Jšœ	œž#˜GJšœ˜Jšœ˜Jšœ˜J˜ J˜J˜J˜J˜J˜J˜J˜J˜Jš œ˜Jš œ˜š œ ˜J˜ J™J™ J™ Jš 4™4J™ J™J™ Jš ¾™¾—J˜ J˜ J˜ š œ œ œ œ œ ˜,J˜ J˜ —Jš ˜˜ Jš D˜DJ˜ J˜ —š J˜JJ˜ Jš L˜LJ˜ J˜ —Jš E˜E˜ Jš 8˜8J˜ J˜ —Jš 4˜4J˜ J˜ Jšœ;˜;J™ Jšœ˜J™šœž˜Jšœ˜Jšœœ˜Jšœ˜—šœž˜"Jšœ˜Jšœœ˜Jšœ˜—šœž>˜DJšœ	œ˜J˜šœ/˜5Jšœ)˜)Jšœ)˜)—š˜Jšœ)˜)Jšœ*˜*—Jšœ1˜1Jšœœœœ˜QJšœ1˜1J™4JšœX˜XJšœ(˜(J™[JšœX˜XJšœ(˜(Jšœ˜—šœžE˜KJšœ˜Jšœœ˜Jšœ˜—šœžA˜GJšœ	œ˜Jšœ™Jšœ˜J˜Jšœ(˜(Jšœ(˜(Jšœ(˜(Jšœ(˜(JšœT˜TJšœœœœ˜QJšœœœ˜5Jšœ]˜]J™5JšœX˜XJšœ(˜(J™[JšœX˜XJšœ(˜(JšœX˜XJšœ(˜(J™[JšœX˜XJšœ(˜(J˜—Jšœœ˜Jšœ˜J˜ —šŸœ	œœœ-˜hJ™Lšœœ˜Lšœž˜ Jšœž˜ Jšœž˜ Jšœž˜ Jšœž˜ Lšœ3ž˜?LšœH™HJšœK™KJšœL™LLšœœœ˜$šœ˜šœœ˜&šœœœ˜(Lšœ0˜0—Lšœ˜Lšœ(˜(L˜—Lšœ˜—Lšœ˜Jšœ˜J˜ —šŸ
œœ)˜9Jšœ8™8Jšœœ˜šœœœ˜#šœœœ˜%šœ)œ˜1Jšœ˜Jšœ(˜(Jšœ˜——Jšœ˜—Jšœ˜Jšœ˜—šŸœœœœ˜GLšœ
œœœ˜Lšœ
œœœ˜Lšœœœ˜Lšœœœœ˜Lšœf˜fLšœ˜L˜ —š
Ÿœœœœœ˜@Lšœ#™#Lšœ
œœ˜Lšœ˜Lšœ˜—šŸœœœœœœœ˜YLšœ&˜&Lšœ˜Lšœ˜L˜—šŸœœ˜Lšœ˜Jšœ-˜-šœ7˜7Lšœ
˜
LšŸœ˜LšŸœ˜Lšœ˜LšÏbœ˜Lšœ"˜"Lšœ)˜)Lšœ˜Lšœ˜Lšœ˜Lšœ˜Lšœ&˜&Lšœ$˜$Lšœ$˜$Lšœ˜Lšœ"˜"—L˜ Lšœ&˜&Lšœ,˜,L˜L˜ —Lšœ˜L˜ Lšœ˜˜ J˜ J˜ —J˜ — …—    ¼  +t  