SVSpheresImpl.mesa
Copyright © 1986 by Xerox Corporation. All rights reserved.
Last edited by Bier on September 23, 1987 10:59:22 pm PDT
Contents: Sphere routines.
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] = {
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.
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;
Transform the line into the natural coordinate system of the sphere.
p ← SVVector3d.Sub[line.base, sphere.center];
d ← line.direction;
[points, hitCount, tangent] ← SphereMeetsLineAux[sphere, p, d, SVUtility.QuadraticFormula];
Transform back to world coordinates.
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;
The sphere is: x^2 + y^2 + z^2 = r^2.
r ← sphere.radius;
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.
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] = {
Use only the positive roots. Preserves the order of the roots, so they will still be in increasing order.
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];
Transform to natural sphere coordinates.
p ← SVVector3d.Sub[p, sphere.center];
[points, hitCount, tangent] ← SphereMeetsLineAux[sphere, p, d, PositiveRoots];
Transform back to world coordinates.
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] = {
If there is a single point of intersection, then tangent will be TRUE iff that point of intersection is a point of tangency.
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) => {
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]
p: Point3d;
normal: Vector3d;
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] = {
We drop a normal from the point onto the sphere and find where it hits.
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.