Skip to content

Instantly share code, notes, and snippets.

Last active Feb 26, 2019
What would you like to do?
Implementations of Ritter's sphere bound algorithm
boundingSphere :: [(Float,Float,Float)] -> ((Float,Float,Float),Float) --takes a list of points in 3D to a pair (center,radius)
boundingSphere points =
case points of --induction on number of points: 0,1,2, or 3+
[] -> ((0,0,0),0)
p:[] -> (p,0)
(p1,p2,p3):(q1,q2,q3):[] -> (((p1+q1)/2,(p2+q2)/2,(p3+q3)/2),dist (p1,p2,p3) (q1,q2,q3))
p:pts -> let
y@(y1,y2,y3) = head $ filter (\pt -> (dist pt p) - 0.1 < (maximum $ map (dist p) pts)) pts
z@(z1,z2,z3) = head $ filter (\pt -> (dist pt y) - 0.1 < (maximum $ map (dist y) pts)) pts
initSphere@(ctr,rad) = (((y1+z1)/2, (y2+z2)/2, (y3+z3)/2), dist y z / 2)
extPts = filter (\pt -> (dist ctr pt) > rad) pts
if length extPts /= 0 then (ctr,maximum $ map (dist ctr) extPts) else initSphere
dist = \(p1,p2,p3) (q1,q2,q3) -> sqrt (p1-q1)^2 + (p2-q2)^2 + (p3-q3)^2
def bounding_sphere(points):
dist = lambda a,b: ((a[0] - b[0])**2 + (a[1] - b[1])**2 + (a[2] - b[2])**2)**0.5 #euclidean metric
x = points[0] #any arbitrary point in the point cloud works
y = max(points,key= lambda p: dist(p,x) ) #choose point y furthest away from x
z = max(points,key= lambda p: dist(p,y) ) #choose point z furthest away from y
bounding_sphere = (((y[0]+z[0])/2,(y[1]+z[1])/2,(y[2]+z[2])/2), dist(y,z)/2) #initial bounding sphere
#while there are points lying outside the bounding sphere, update the sphere by growing it to fit
exterior_points = [p for p in points if dist(p,bounding_sphere[0]) > bounding_sphere[1] ]
while ( len(exterior_points) > 0 ):
pt = exterior_points.pop()
if (dist(pt, bounding_sphere[0]) > bounding_sphere[1]):
bounding_sphere = (bounding_sphere[0],dist(pt,bounding_sphere[0])) #grow the sphere to capture new point
return bounding_sphere

This comment has been minimized.

Copy link
Owner Author

@paultsw paultsw commented May 1, 2014

N.B. Ritter's bounding algorithm (see is not optimal, i.e. it does not return the sphere of smallest radius that bounds all points; in general it is anywhere from five to twenty percent off from the minimal radius sphere. Further, the implementations are slightly different due to the addition of a small additive factor in the Haskell implementation for the sake of "leeway" when filtering out the list of points; this is to make sure that certain floats "close enough" to the boundary don't get discounted.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment