<> <> <> 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] = { <> << 2 2 2 2 2 2 F(x,y,z) = y z + x z - 2 x y z + x y = 0 >> <> 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; <> <<>> <<>> << 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 <> 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 <<[realRootsOneOrigin, rootCount] _ Roots.RealQuartic[a4, a3, a2, a1, a0];>> <> <> 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.