File: SVBoundSphereImpl.mesa
Last edited by: Eric Bier on August 1, 1985 1:26:54 am PDT
Copyright © 1984 by Xerox Corporation. All rights reserved.
Contents: Computes a sphere which is guaranteed to bound a given polyhedron. Implements operations on such spheres.
DIRECTORY
CoordSys,
CSGGraphics,
Matrix3d,
RealFns,
Imager,
SV3d,
SVBoundBox,
SVBoundSphere,
SVModelTypes,
SVPolygon3d,
SVRayTypes,
SVVector3d;
SVBoundSphereImpl:
CEDAR
PROGRAM
IMPORTS CoordSys, CSGGraphics, Matrix3d, RealFns, SVBoundBox, SVPolygon3d, SVVector3d
EXPORTS SVBoundSphere =
BEGIN
BoundHedron: TYPE = SVModelTypes.BoundHedron;
Camera: TYPE = SVModelTypes.Camera;
CoordSystem: TYPE = SVModelTypes.CoordSystem;
Matrix4by4: TYPE = SV3d.Matrix4by4;
Point3d: TYPE = SV3d.Point3d;
Poly3d: TYPE = SV3d.Poly3d;
Ray: TYPE = SVRayTypes.Ray;
Vector: TYPE = SV3d.Vector;
BoundSphere: TYPE = REF BoundSphereObj;
BoundSphereObj: TYPE = SVModelTypes.BoundSphereObj;
gCirclePolygon: Poly3d;
PI: REAL = 3.14159265358979;
Init:
PROC [] = {
gCirclePolygon ← PolyOnUnitCircle[20];
};
PolyOnUnitCircle:
PROC [edges:
NAT]
RETURNS [poly: Poly3d] ~
TRUSTED {
Generate a regular polygon whose vertices sit on the unit circle, centered on the origin, in the z=0 plane.
theta: REAL ← 0.0;
realEdges: REAL ← edges;
delta: REAL ← 2.0*PI/realEdges;
x, y: REAL;
poly ← SVPolygon3d.CreatePoly[edges];
FOR i:
NAT
IN [0..edges)
DO
x ← RealFns.Sin[theta];
y ← RealFns.Cos[theta];
poly ← SVPolygon3d.AddPolyPoint[poly, [x, y, 0]];
theta ← theta + delta;
ENDLOOP;
};
CopyBoundSphere:
PUBLIC
PROC [boundSphere: BoundSphere]
RETURNS [newBS: BoundSphere] = {
newBS ←
NEW[BoundSphereObj ← [
boundSphere.center,
boundSphere.radius,
boundSphere.radiusSquared]];
};
BoundSphereFromBoundHedron:
PUBLIC
PROC [bh: BoundHedron, worldCS: CoordSystem, localCS: CoordSystem]
RETURNS [boundSphere: BoundSphere] =
TRUSTED {
For now we just use a simple heuristic. Use the center of mass as the center of the sphere. Use the distance from this center to the farthest point as the radius. Remember to convert all coordinates to WORLD.
cmWORLD, bhPtWORLD, localCM: Point3d;
maxIndex: NAT;
maxDistSquared, distSquared, radius: REAL;
IF bh = NIL THEN {boundSphere ← NIL; RETURN};
localCM ← SVBoundBox.CenterOfMassBoundHedron[bh];
cmWORLD ← Matrix3d.Update[localCS.wrtWorld, localCM];
maxIndex ← 0;
bhPtWORLD ← Matrix3d.Update[localCS.wrtWorld, bh[0]];
maxDistSquared ← SVVector3d.MagnitudeSquared[SVVector3d.Sub[cmWORLD, bhPtWORLD]];
FOR i:
NAT
IN[1..bh.len)
DO
bhPtWORLD ← Matrix3d.Update[localCS.wrtWorld, bh[i]];
distSquared ← SVVector3d.MagnitudeSquared[SVVector3d.Sub[cmWORLD, bhPtWORLD]];
IF distSquared > maxDistSquared
THEN {
maxIndex ← i;
maxDistSquared ← distSquared;
};
ENDLOOP;
radius ← RealFns.SqRt[maxDistSquared];
boundSphere ←
NEW[BoundSphereObj ←
[center: cmWORLD, radius: radius, radiusSquared: radius*radius]];
};
BoundSphereFromValues:
PUBLIC
PROC [center
WORLD: Point3d, radius:
REAL]
RETURNS [boundSphere: BoundSphere] = {
boundSphere ← NEW[BoundSphereObj ← [centerWORLD, radius, radius*radius]];
};
UnionCombineBoundSpheres:
PUBLIC
PROC [bs1, bs2: BoundSphere]
RETURNS [newBS: BoundSphere] =
TRUSTED {
A sphere which contains bs1 and bs2 will have diameter r1 + ||c2-c1|| +r2 (or 2r1 if bs1 completely contains bs2 and vice-versa). One diameter will contain the points c1-r1d and c2+r2d where d is the unit vector in the c2-c1 direction. Hence, the center is at (c1-r1d + c2+r2d)/2 = (c1 + c2 + (r2-r1)d)/2.
d, unitD: Vector;
almostZero: REAL ← 1.0e-12;
magD: REAL;
IF bs1 = NIL OR bs2 = NIL THEN {newBS ← NIL; RETURN};
d ← SVVector3d.Sub[bs2.center, bs1.center];
magD ← SVVector3d.Magnitude[d];
If one sphere contains the other, handle as a special case and RETURN
IF bs1.radius > bs2.radius
THEN {
IF bs1.radius - magD >= bs2.radius THEN {newBS ← CopyBoundSphere[bs1]; RETURN};
}
ELSE {
IF bs2.radius - magD >= bs1.radius THEN {newBS ← CopyBoundSphere[bs2]; RETURN};
};
There is a problem if the two spheres are almost identical.
IF
ABS[magD] < almostZero
THEN {
IF bs1.radius > bs2.radius THEN {newBS ← CopyBoundSphere[bs1]; RETURN}
ELSE {newBS ← CopyBoundSphere[bs2]; RETURN}
};
Normal case.
newBS ← NEW[BoundSphereObj];
unitD ← SVVector3d.Scale[d, 1.0/magD];
newBS.radius ← (bs1.radius + magD + bs2.radius)/2.0;
newBS.radiusSquared ← newBS.radius*newBS.radius;
newBS.center ← SVVector3d.Add[bs1.center, bs2.center];
newBS.center ← SVVector3d.Add[newBS.center, SVVector3d.Scale[unitD, bs2.radius - bs1.radius]];
newBS.center ← SVVector3d.Scale[newBS.center, 0.5];
};
IntersectionCombineBoundSpheres:
PUBLIC
PROC [bs1, bs2: BoundSphere]
RETURNS [newBS: BoundSphere] =
TRUSTED {
Pretend that bs1 is centered on [0,0] and bs2 on [c,0] then they intersect, if at all, at the points [x, +-SqRt[r2 - x2]], where x = (r2-R2+c2)/2c. Hence, for general bs1 and bs2, their intersection is bounded by the sphere centered at (x/c)(bs1.center-bs2.center) + bs1.center, with radius SqRt[r2 - x2].
c, cSquared, x, y: REAL;
o2minuso1: Vector;
IF bs1 = NIL AND bs2 = NIL THEN {newBS ← NIL; RETURN};
IF bs1 = NIL THEN {newBS ← CopyBoundSphere[bs2]; RETURN};
IF bs2 = NIL THEN {newBS ← CopyBoundSphere[bs1]; RETURN};
o2minuso1 ← SVVector3d.Sub[bs2.center, bs1.center];
cSquared ← SVVector3d.MagnitudeSquared[o2minuso1];
c ← RealFns.SqRt[cSquared];
IF c > (bs1.radius + bs2.radius)
THEN { -- No intersection. Return a radius zero sphere.
newBS ← NEW[BoundSphereObj];
newBS.center ← bs1.center;
newBS.radius ← 0.0;
newBS.radiusSquared ← 0.0;
RETURN;
};
IF (c + bs1.radius) < bs2.radius THEN { -- bs1 is contained in bs2
newBS ← CopyBoundSphere[bs1]; RETURN;
};
IF (c + bs2.radius) < bs1.radius THEN { -- bs2 is contained in bs1
newBS ← CopyBoundSphere[bs2];
};
newBS ← NEW[BoundSphereObj];
x ← (bs1.radiusSquared - bs2.radiusSquared + cSquared)/(2.0*c);
y ← RealFns.SqRt[bs1.radiusSquared - x*x];
newBS.center ← SVVector3d.Add[SVVector3d.Scale[o2minuso1, x/c], bs1.center];
newBS.radius ← y;
newBS.radiusSquared ← y*y;
};
DifferenceCombineBoundSpheres:
PUBLIC
PROC [bs1, bs2: BoundSphere]
RETURNS [newBS: BoundSphere] =
TRUSTED {
For most cases, just use the first (positive) sphere. It is hard to do any better than this.
IF bs1 = NIL THEN {newBS ← NIL; RETURN};
newBS ← CopyBoundSphere[bs1];
};
DrawBoundSphere:
PUBLIC
PROC [dc: Imager.Context, boundSphere: BoundSphere, camera: Camera] =
TRUSTED {
Find a coordinate system S whose z-axis points from sphere center to eyepoint, origin on the sphere center and zy plane containing the y-axis of CAMERA. We describe S in terms of CAMERA coordinates. Scale SCAMERA by the radius of the boundSphere and draw absolutely (since all points are already in CAMERA coords).
circleMat, unitCircleMat: Matrix4by4;
circlePoly: Poly3d;
rad: REAL ← boundSphere.radius;
centerToEyepoint: Vector;
centerCAMERA: Point3d;
worldCS: CoordSystem;
eyepointCAMERA: Point3d ← [0, 0, camera.focalLength];
worldCS ← camera.coordSys.parent;
centerCAMERA ← CoordSys.FromCSToCS[boundSphere.center, worldCS, camera.coordSys];
centerToEyepoint ← SVVector3d.VectorFromPoints[centerCAMERA, eyepointCAMERA];
unitCircleMat ← MakeAlignedMat[centerToEyepoint, centerCAMERA];
circleMat ← Matrix3d.LocalScale[unitCircleMat, rad, rad, rad];
circlePoly ← SVPolygon3d.TransformByMat[gCirclePolygon, circleMat];
CSGGraphics.DrawPolygonAbsolute[dc, circlePoly, 1, butt, camera];
};
MakeAlignedMat:
PRIVATE
PROC [zaxis: Vector, origin: Point3d]
RETURNS [mat: Matrix4by4] =
TRUSTED {
Create a Matrix4by4 with origin at origin whose z axis is parallel to zaxis and whose x axis is orthogonal to both zaxis and the y axis of the reference coordinate frame (CAMERA).
xAxis: Vector;
yAxisOfCamera: Vector ← [0,1,0];
IF SVVector3d.Parallel[yAxisOfCamera, zaxis] THEN xAxis ← [1,0,0]
ELSE xAxis ← SVVector3d.CrossProduct[yAxisOfCamera, zaxis];
mat ← Matrix3d.MakeMatFromZandXAxis[zaxis, xAxis, origin];
};
Init[];
END.