[Indigo]<Solidviews>SolidModeler>SVBoundSphereImpl.mesa!7

File: SVBoundSphereImpl.mesa

Last edited by: Eric Bier on May 4, 1987 11:36:29 pm PDT

Contents: Computes a sphere which is guaranteed to bound a given polyhedron. Implements operations on such spheres.

DIRECTORY

CoordSys, SVGraphics, Matrix3d, RealFns, Imager, SV3d, SVBasicTypes, SVBoundBox, SVBoundSphere, SVModelTypes, SVPolygon3d, SVRayTypes, SVVector3d;

SVBoundSphereImpl: CEDAR PROGRAM

IMPORTS CoordSys, SVGraphics, Matrix3d, RealFns, SVBoundBox, SVPolygon3d, SVVector3d

EXPORTS SVBoundSphere =

BEGIN

BoundHedron: TYPE = SVBasicTypes.BoundHedron;

BoundSphere: TYPE = REF BoundSphereObj;

BoundSphereObj: TYPE = SVBasicTypes.BoundSphereObj;

Camera: TYPE = SVModelTypes.Camera;

CoordSystem: TYPE = SVModelTypes.CoordSystem;

Matrix4by4: TYPE = SV3d.Matrix4by4;

Point3d: TYPE = SV3d.Point3d;

Poly3d: TYPE = SV3d.Poly3d;

Ray: TYPE = SVRayTypes.Ray;

Vector3d: TYPE = SV3d.Vector3d;

gCirclePolygon: Poly3d;

PI: REAL = 3.14159265358979;

Init: PROC [] = {

gCirclePolygon ← PolyOnUnitCircle[20];

};

PolyOnUnitCircle: PROC [edges: NAT] RETURNS [poly: Poly3d] ~ {

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] = {

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, cmlocal: Point3d;

maxIndex: NAT;

maxDistSquared, distSquared, radius: REAL;

localWorld: Matrix4by4 ← CoordSys.WRTWorld[localCS];

IF bh = NIL THEN {boundSphere ← NIL; RETURN};

cmlocal ← SVBoundBox.CenterOfMassBoundHedron[bh];

cmWORLD ← Matrix3d.Update[cmlocal, localWorld];

maxIndex ← 0;

bhPtWORLD ← Matrix3d.Update[bh[0], localWorld];

maxDistSquared ← SVVector3d.MagnitudeSquared[SVVector3d.Sub[cmWORLD, bhPtWORLD]];

FOR i: NAT IN[1..bh.len) DO

bhPtWORLD ← Matrix3d.Update[bh[i], localWorld];

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 [centerWORLD: Point3d, radius: REAL] RETURNS [boundSphere: BoundSphere] = {

boundSphere ← NEW[BoundSphereObj ← [centerWORLD, radius, radius*radius]];

};

UnionCombineBoundSpheres: PUBLIC PROC [bs1, bs2: BoundSphere] RETURNS [newBS: BoundSphere] = {

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: Vector3d;

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] = {

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: Vector3d;

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] = {

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] = {

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;

centerToEyepoint: Vector3d;

centerCAMERA: Point3d;

worldCS: CoordSystem;

eyepointCAMERA: Point3d ← [0, 0, camera.focalLength];

IF boundSphere = NIL THEN RETURN;

rad ← boundSphere.radius;

worldCS ← CoordSys.Parent[camera.coordSys];

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];

SVGraphics.DrawPolygonAbsolute[dc, circlePoly, 1, butt, camera];

};

MakeAlignedMat: PRIVATE PROC [zaxis: Vector3d, origin: Point3d] RETURNS [mat: Matrix4by4] = {

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: Vector3d;

yAxisOfCamera: Vector3d ← [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.