Instantly share code, notes, and snippets.

# paultsw/boundingSphere.hs

Last active February 26, 2019 02:29
Show Gist options
• Save paultsw/074b612d4e5980a18a94 to your computer and use it in GitHub Desktop.
Implementations of Ritter's sphere bound algorithm
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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

### 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.