{{ message }}

Instantly share code, notes, and snippets.

# paultsw/boundingSphere.hs

Last active Feb 26, 2019
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 in if length extPts /= 0 then (ctr,maximum \$ map (dist ctr) extPts) else initSphere where 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 - b)**2 + (a - b)**2 + (a - b)**2)**0.5 #euclidean metric x = points #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+z)/2,(y+z)/2,(y+z)/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) > bounding_sphere ] while ( len(exterior_points) > 0 ): pt = exterior_points.pop() if (dist(pt, bounding_sphere) > bounding_sphere): bounding_sphere = (bounding_sphere,dist(pt,bounding_sphere)) #grow the sphere to capture new point return bounding_sphere

### paultsw commented May 1, 2014

 N.B. Ritter's bounding algorithm (see https://en.wikipedia.org/wiki/Bounding_sphere#Ritter.27s_bounding_sphere) 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.