Implementations of Ritter's sphere bound algorithm

boundingSphere.hs

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

boundingSphere.py

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

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.