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.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] = { 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; 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]] = { 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]; 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) => { 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. dSVSpheresImpl.mesa Copyright c 1986 by Xerox Corporation. All rights reserved. Last edited by Bier on September 23, 1987 10:59:22 pm PDT Contents: Sphere routines. "sphere" now has center at "point" and radius of "radius". "sphere" now passes through the four points. If the points are roughly coplanar, creates a large sphere passing through p0, p1 and p2. Returns a sphere which passes through the three points. If the points are roughly collinear, creates a large sphere passing through p0, p1 and p2. Transform the line into the natural coordinate system of the sphere. Transform back to world coordinates. The sphere is: x^2 + y^2 + z^2 = r^2. The line has the form: x = p1 + t*d1 y = p2 + t*d2 z = p3 + t*d3 for -# < t < #. hence, at the intersection of ray and sphere, t satisfies: (p1 + t*d1)^2 + (p2 + t*d2)^2 + (p3 + t*d3)^2 = r^2 which reduces to: t^2 * d.d + t * 2 * d.p + p.p - r^2 = 0. Use only the positive roots. Preserves the order of the roots, so they will still be in increasing order. Transform to natural sphere coordinates. Transform back to world coordinates. If there is a single point of intersection, then tangent will be TRUE iff that point of intersection is a point of tangency. "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. [Artwork node; type 'ArtworkInterpress on' to command tool] The two spheres overlap. We expect two roots. In the picture below, point C is on the circle of intersection. Segment AB is the segment O1O2 in the picture above. We wish to find the altitude h of the triangle. b = r1 and a = r2. From Heron's formula for the area of a triangle, area K = sqrt(s(s-a)(s-b)(s-c)), where s = (a+b+ c)/2. We also know that the area K = 0.5ch. So h = 2K/c. This gives us the radius of the intersection circle. The x coordinate is the point where the altitude hits the base, which we get from the Pythagorean Theorem. [Artwork node; type 'ArtworkInterpress on' to command tool] normal: Vector3d; We drop a normal from the point onto the sphere and find where it hits. ΚG˜J˜Icodešœ™Kšœ Οmœ1™Kšœ žœ ˜K˜—šŸ œžœžœ˜6Kšœ ˜ K˜—K˜šŸœžœžœžœ˜ZKšœ:™:Kšœ˜Kšœ˜Kšœ%˜%K˜—š Ÿœžœžœžœžœ˜^Kšœ žœ*˜6K˜—K˜š Ÿœžœžœ,žœ žœ˜hK™‡Kšœ žœ˜Kšžœ<˜BK˜—š Ÿœžœžœžœžœ˜_Kšœ“™“Kšœ žœ˜Kšžœ8˜>K˜—K˜šŸœžœžœ!žœ žœžœ%žœ˜‹Kšœ ˜ Kšœ ˜ KšΟbD™DKšœ-˜-Kšœ˜Kšœ[˜[Kš $™$šžœžœžœž˜Kšœ5˜5Kšžœ˜—K˜K˜—Kšœžœžœ žœžœ žœžœžœ žœ˜aK˜š ŸœžœFŸžœ žœžœ%žœ˜­Kšœžœ˜Iprocšœžœžœžœ˜Lšœž˜ Lšœž˜ Lšœž˜ Lšœž˜ Lšœžœ˜ Lš %™%Lšœ˜Lš ™Lš ™Lš ™Lš ™Lš Πbm ‘ <™JLš 3™3Lš ™Lš (™(L˜Lšœ!˜!Lšœ ˜ Lšœ ˜ Lšœ ˜ L˜Lšœ0˜0L™šžœ ž˜šœ˜Kšœ+žœ˜1Lšœ˜—šœ˜Lšœ˜Lšœ˜Lšœ˜Lšœ˜Lšœ žœ˜Lšœ˜—šœ˜Lšœ ˜ Lšœ˜Lšœ˜Lšœ˜Lšœ˜Lšœ˜Lšœ˜Lšœ žœ˜Lšœ˜—Lšžœžœ˜—K˜K˜—šŸ œžœ žœžœ žœžœžœ žœ˜\J™jLšœžœ˜ Lšœ9˜9Lšžœžœžœ˜šžœž˜šžœžœ˜šžœžœžœž˜!Lšœ˜—Lšžœ˜Lšœ˜L˜—Lšžœ ˜—Lšžœ˜Jšœ˜—šŸœžœžœžœ žœžœžœžœ&žœ˜©K˜ K˜ Kšœ ˜ Kš (™(Kšœ%˜%KšœN˜NKš $™$šžœžœžœž˜Kšœ˜Kšœ5˜5Kšžœ˜—K˜K˜—šŸœžœžœ!žœ žœžœ%žœ˜‹KšœAžœ7™|Kšœ žœžœ ˜$Kšœ˜KšœF˜FKšœ ˜ šžœžœžœž˜šžœ0žœ˜8Kšœ˜Kšœ!˜!K˜—Kšžœ˜—K˜K˜—šŸ œžœžœžœ˜6Kšœ žœ ˜Kš žœžœžœžœžœžœ˜gK˜K˜—šŸœžœžœžœ+˜mKšœ•™•Icaption˜Mšœ žœ˜*Kšœ'žœ˜,šžœ-žœΟc˜KKšžœ!žœžœžœ˜6Kšžœžœžœ˜K˜—Kšœ8˜8Kšœ)˜)Kšœ/˜/Kšœžœ"˜4Icenter–G49.03611 mm topLeading 49.03611 mm topIndent 1.411111 mm bottomLeading •Bounds:0.0 mm xmin 0.0 mm ymin 67.73333 mm xmax 46.21388 mm ymax •Artwork Interpress• InterpressΩInterpress/Xerox/3.0  f j k j‘₯“ΔWB € ¨  `#‘£§G ’ ¨Δ™ΔώS ’ ¨‘‘¨ΔWB € ¨ r j‘₯“ΔWB € ¨‘‘¨‘―“‘‘¨’°“ ‘¨ή™Hή‘“‘‘‘‘™‘‘¨ή™Hή‘“‘Έ‘―“‘―“‘‘¨’°“ ‘¨H™€H‘“‘‘‘‘™‘‘¨H™€H‘“‘Έ‘―“‘―“‘‘¨’°“ΔšHΔ‘Ή‘―“‘―“‘‘¨’°“Δs Δa‘Ίd‘Ή‘―“ΕXeroxΕ PressFontsΕ Helvetica-mrr£‘ “ͺ € ” •  —Δ:ƒ/ŠΑr1–  —e ŠΑr2–  —ŠΑO1–  —k ŠΑO2–‘―“‘‘¨’°“d‘Ή‘―“ k ι k gšœ=™=Kšžœ ž˜Kšœ$žœ’˜Nšžœ5’˜KKšœ˜Kšœ žœ˜Kšœ4˜4Kšœ\˜\Kšœ˜KšœM˜MK˜—šžœ4˜6Kšœyžœžœ›™ͺN–G81.13889 mm topLeading 81.13889 mm topIndent 1.411111 mm bottomLeading –:0.0 mm xmin 0.0 mm ymin 112.8889 mm xmax 78.31667 mm ymax – Interpress–ŽInterpress/Xerox/3.0  f j k j‘₯“ΔWB € ¨  ΰ~‘£ί» ’ ¨  ‘ ’ ¨‘‘¨ΔWB € ¨ r j‘₯“ΔWB € ¨‘‘¨‘―“‘‘¨’°“x›πL‘Ή‘―“‘‘¨’°“πL…›‘Ή‘―“‘‘¨’°“…›x›‘Ή‘―“ΕXeroxΕ PressFontsΕ Helvetica-mrr£‘ “ͺ € ” •  —@ώŠΑa–  —©ϊŠΑb–  —Δ>_0ΔΣ?γŠΑc–‘―“‘‘¨’°“π›πL‘Ή‘―“  —ηĊΑh–  —Έ©ŠΑm–  —rŠΑA–  —Δ0G”ŠΑB–  —ΔQί?SŠΑC– k ι k gšœ=™=K˜ K™Kšœ žœ˜Kšœžœ˜Kšœžœ˜K˜Kšœ4˜4Kšœ˜Kšœ<˜<šžœžœ˜Kšœ˜KšœC˜CK˜—šžœ˜Kšœ˜KšœC˜CK˜—Kšœ žœ˜K˜K˜KšœA˜AK˜—šžœ4˜6J˜Kšœ žœ˜šžœ!žœ˜)Kšœ4˜4Kšœ\˜\K˜—šžœ˜Kšœ5˜5Kšœ\˜\K˜—K˜KšœM˜MK˜—šœ˜Kšœ˜K˜Kšœžœ˜K˜—Kšžœžœ˜Kšœ˜K˜—š Ÿœžœžœ žœžœ˜VKšžœ9˜?K˜—š Ÿœžœžœžœžœ˜OKšœžœ#˜*K˜—šŸœžœžœ žœ˜gKšœG™GKšœ žœ ˜K˜Kšœžœ˜Kšœ2˜2Kšœ7˜7šžœžœ˜$KšœA˜AKšžœ˜K˜—Kšœm˜mK˜K˜—K˜Kšžœ˜K˜J˜—…—v4!