<> <> <> <> <<>> DIRECTORY Feedback, RealFns, Rope, SV3d, SVLines3d, SVPolygon3d, SVRay, SVSceneTypes, SVSpheres, SVUtility, SVVector3d; SVSpheresImpl: CEDAR PROGRAM IMPORTS Feedback, RealFns, SVLines3d, SVPolygon3d, SVRay, SVUtility, SVVector3d EXPORTS SVSpheres = BEGIN Circle3d: TYPE = REF Circle3dObj; Circle3dObj: TYPE = SV3d.Circle3dObj; Edge3d: TYPE = SV3d.Edge3d; Line3d: TYPE = SV3d.Line3d; Point3d: TYPE = SV3d.Point3d; Ray: TYPE = SVSceneTypes.Ray; Sphere: TYPE = REF SphereObj; SphereObj: TYPE = SV3d.SphereObj; Vector3d: TYPE = SV3d.Vector3d; Problem: SIGNAL [msg: Rope.ROPE] = Feedback.Problem; CreateEmptySphere: PUBLIC PROC [] RETURNS [sphere: Sphere] = { sphere _ NEW[SphereObj]; }; CopySphere: PUBLIC PROC [from: Sphere, to: Sphere] = { to^ _ from^; }; FillSphereFromPointAndRadius: PUBLIC PROC [pt: Point3d, radius: REAL, sphere: Sphere] = { <<"sphere" now has center at "point" and radius of "radius".>> sphere.center _ pt; sphere.radius _ radius; sphere.radiusSquared _ radius*radius; }; SphereFromPointAndRadius: PUBLIC PROC [pt: Point3d, radius: REAL] RETURNS [sphere: Sphere] = { sphere _ NEW[SphereObj _ [pt, radius, radius*radius]]; }; FillSphereFrom4Points: PUBLIC PROC [p0, p1, p2, p3: Point3d, sphere: Sphere] RETURNS [planar: BOOL] = { <<"sphere" now passes through the four points. If the points are roughly coplanar, creates a large sphere passing through p0, p1 and p2.>> planar _ FALSE; SIGNAL Problem[msg: "Not yet implemented: FillSphereFrom4Points"]; }; SphereFrom4Points: PUBLIC PROC [p0, p1, p2: Point3d] RETURNS [sphere: Sphere, planar: BOOL] = { <> planar _ FALSE; SIGNAL Problem[msg: "Not yet implemented: SphereFrom4Points"]; }; SphereMeetsLine: PUBLIC PROC [sphere: Sphere, line: Line3d] RETURNS [points: ARRAY [1..2] OF Point3d, hitCount: [0..2], tangent: BOOL] = { p: Point3d; d: Vector3d; <> p _ SVVector3d.Sub[line.base, sphere.center]; d _ line.direction; [points, hitCount, tangent] _ SphereMeetsLineAux[sphere, p, d, SVUtility.QuadraticFormula]; <> FOR i: NAT IN [1..hitCount] DO points[i] _ SVVector3d.Add[points[i], sphere.center]; ENDLOOP; }; RootFinderProc: TYPE = PROC [a, b, c: REAL] RETURNS [roots: ARRAY[1..2] OF REAL, rootCount: NAT]; SphereMeetsLineAux: PROC [sphere: Sphere, p: Point3d, d: Vector3d, rootFinder: RootFinderProc] RETURNS [points: ARRAY [1..2] OF Point3d, hitCount: [0..2], tangent: BOOL] = { r: REAL; tArray: ARRAY[1..2] OF REAL; dp: REAL; dp2: REAL; dd: REAL; pp: REAL; t1, t2: REAL; <> r _ sphere.radius; <> << x = p1 + t*d1>> << y = p2 + t*d2>> << z = p3 + t*d3>> <> <<(p1 + t*d1)^2 + (p2 + t*d2)^2 + (p3 + t*d3)^2 = r^2>> <> <> dp _ SVVector3d.DotProduct[d,p]; dp2 _ dp*dp; dd _ SVVector3d.DotProduct[d,d]; pp _ SVVector3d.DotProduct[p,p]; [tArray, hitCount] _ rootFinder[dd,2*dp,pp-r*r]; <<>> SELECT hitCount FROM 0 => { points[1] _ points[2] _ [0,0,0]; tangent _ FALSE; }; 1 => { t1 _ tArray[1]; points[1][1] _ p[1] + d[1]*t1; points[1][2] _ p[2] + d[2]*t1; points[1][3] _ p[3] + d[3]*t1; tangent _ TRUE; }; 2 => { t1 _ tArray[1]; t2 _ tArray[2]; points[1][1] _ p[1] + d[1]*t1; points[1][2] _ p[2] + d[2]*t1; points[1][3] _ p[3] + d[3]*t1; points[2][1] _ p[1] + d[1]*t2; points[2][2] _ p[2] + d[2]*t2; points[2][3] _ p[3] + d[3]*t2; tangent _ FALSE; }; ENDCASE => ERROR; }; PositiveRoots: PROC [a, b, c: REAL] RETURNS [roots: ARRAY[1..2] OF REAL, rootCount: NAT] = { <> i: NAT _ 1; [roots, rootCount] _ SVUtility.QuadraticFormula[a, b, c]; IF rootCount = 0 THEN RETURN; WHILE i <= rootCount DO IF roots[i] <= 0 THEN { FOR j: NAT IN [i..rootCount-1] DO roots[j] _ roots[j+1]; ENDLOOP; rootCount _ rootCount - 1; } ELSE {i _ i + 1}; ENDLOOP; }; SphereMeetsRay: PUBLIC PROC [sphere: Sphere, ray: Ray] RETURNS [points: ARRAY [1..2] OF Point3d, normals: ARRAY [1..2] OF Vector3d, hitCount: [0..2], tangent: BOOL] = { p: Point3d; d: Vector3d; [p, d] _ SVRay.GetWorldRay[ray]; <> p _ SVVector3d.Sub[p, sphere.center]; [points, hitCount, tangent] _ SphereMeetsLineAux[sphere, p, d, PositiveRoots]; <> FOR i: NAT IN [1..hitCount] DO normals[i] _ points[i]; points[i] _ SVVector3d.Add[points[i], sphere.center]; ENDLOOP; }; SphereMeetsEdge: PUBLIC PROC [sphere: Sphere, edge: Edge3d] RETURNS [points: ARRAY [1..2] OF Point3d, hitCount: [0..2], tangent: BOOL] = { <> testPoints: ARRAY [1..2] OF Point3d; testCount: [0..2]; [testPoints, testCount, tangent] _ SphereMeetsLine[sphere, edge.line]; hitCount _ 0; FOR i: NAT IN [1..testCount] DO IF SVLines3d.LinePointOnEdge[testPoints[i], edge] THEN { hitCount _ hitCount + 1; points[hitCount] _ testPoints[i]; }; ENDLOOP; }; RatherClose: PROC [p1, p2: Point3d] RETURNS [BOOL] = { epsilon: REAL = 1.0e-5; RETURN[ABS[p1[1] - p2[1]] < epsilon AND ABS[p1[2] - p2[2]] < epsilon AND ABS[p1[3] - p2[3]] < epsilon]; }; SphereMeetsSphere: PUBLIC PROC [sphere1, sphere2: Sphere] RETURNS [circle: Circle3d, dimension: [-1..2]] = { <<"dimension" is the dimension of the intersection: -1 for a miss, 0 for a point of tangency, 1 for a circle of intersection, 2 for coincident spheres.>> o1ToO2, o1ToO2Hat: Vector3d; epsilon: REAL = SVUtility.epsilonInPoints; magO1ToO2, outerTangent, innerTangent: REAL; IF RatherClose[sphere1.center, sphere2.center] THEN { -- concentric spheres IF sphere1.radius = sphere2.radius THEN RETURN[NIL, 2] ELSE RETURN [NIL, -1]; }; o1ToO2 _ SVVector3d.Sub[sphere2.center, sphere1.center]; magO1ToO2 _ SVVector3d.Magnitude[o1ToO2]; outerTangent _ sphere1.radius + sphere2.radius; innerTangent _ ABS[sphere1.radius - sphere2.radius]; << [Artwork node; type 'ArtworkInterpress on' to command tool] >> SELECT magO1ToO2 FROM > outerTangent+epsilon => {circle _ NIL; dimension _ -1}; -- spheres far apart IN [outerTangent-epsilon .. outerTangent+epsilon] => { -- point of tangency dimension _ 0; circle _ NEW[Circle3dObj]; o1ToO2Hat _ SVVector3d.Scale[o1ToO2, 1.0/magO1ToO2]; circle.origin _ SVVector3d.Add[sphere1.center, SVVector3d.Scale[o1ToO2Hat, sphere1.radius]]; circle.radius _ 0.0; circle.plane _ SVPolygon3d.PlaneFromPointAndNormal[circle.origin, o1ToO2Hat]; }; IN (innerTangent+epsilon .. outerTangent-epsilon) => { <> << [Artwork node; type 'ArtworkInterpress on' to command tool] >> p: Point3d; <> s, h, m: REAL; b: REAL _ sphere1.radius; a: REAL _ sphere2.radius; dimension _ 1; o1ToO2Hat _ SVVector3d.Scale[o1ToO2, 1.0/magO1ToO2]; s _ 0.5*(b + a + magO1ToO2); h _ 2.0*RealFns.SqRt[s*(s-magO1ToO2)*(s-b)*(s-a)]/magO1ToO2; IF b > a THEN { m _ RealFns.SqRt[b*b - h*h]; p _ SVVector3d.Add[sphere1.center, SVVector3d.Scale[o1ToO2Hat, m]]; } ELSE { m _ RealFns.SqRt[a*a - h*h]; p _ SVVector3d.Sub[sphere2.center, SVVector3d.Scale[o1ToO2Hat, m]]; }; circle _ NEW[Circle3dObj]; circle.origin _ p; circle.radius _ h; circle.plane _ SVPolygon3d.PlaneFromPointAndNormal[p, o1ToO2Hat]; }; IN [innerTangent-epsilon .. innerTangent+epsilon] => { dimension _ 0; circle _ NEW[Circle3dObj]; IF sphere1.radius > sphere2.radius THEN { o1ToO2Hat _ SVVector3d.Scale[o1ToO2, 1.0/magO1ToO2]; circle.origin _ SVVector3d.Add[sphere1.center, SVVector3d.Scale[o1ToO2Hat, sphere1.radius]]; } ELSE { o1ToO2Hat _ SVVector3d.Scale[o1ToO2, -1.0/magO1ToO2]; circle.origin _ SVVector3d.Add[sphere2.center, SVVector3d.Scale[o1ToO2Hat, sphere2.radius]]; }; circle.radius _ 0.0; circle.plane _ SVPolygon3d.PlaneFromPointAndNormal[circle.origin, o1ToO2Hat]; }; < innerTangent-epsilon => { dimension _ -1; circle.radius _ 0.0; circle.plane _ NIL; }; ENDCASE => ERROR; }; SignedSphereDistance: PUBLIC PROC [pt: Point3d, sphere: Sphere] RETURNS [d: REAL] = { RETURN[SVVector3d.Distance[pt, sphere.center] - sphere.radius]; }; SphereDistance: PUBLIC PROC [pt: Point3d, sphere: Sphere] RETURNS [d: REAL] = { d _ ABS[SignedSphereDistance[pt, sphere]]; }; PointProjectedOntoSphere: PUBLIC PROC [pt: Point3d, sphere: Sphere] RETURNS [projectedPt: Point3d] = { <> epsilon: REAL = 1.0e-5; originToPoint: Vector3d; magOriginToPoint: REAL; originToPoint _ SVVector3d.Sub[pt, sphere.center]; magOriginToPoint _ SVVector3d.Magnitude[originToPoint]; IF magOriginToPoint < epsilon THEN { projectedPt _ SVVector3d.Add[sphere.center, [sphere.radius,0,0]]; RETURN; }; projectedPt _ SVVector3d.Add[sphere.center, SVVector3d.Scale[originToPoint, sphere.radius/magOriginToPoint]]; }; END.