Ritter's algorithm

ritter.py

from scipy.stats import uniform
import numpy as np

def dist(a, b) -> float:
    return  ((a[0] - b[0])**2 + (a[1] - b[1])**2 + (a[2] - b[2])**2)**0.5
    
N = 10000


def bounding_sphere(points):
    """
    From https://gist.github.com/paultsw/074b612d4e5980a18a94
    """
    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
   
def bounding_sphere_ritter(points):
    # Initial guess, same as before
    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
    center, radius = (((y[0]+z[0])/2,(y[1]+z[1])/2,(y[2]+z[2])/2), dist(y,z)/2)    #initial bounding sphere
    
    # Note this doesn't use the radius^2 optimization that Ritter uses
    for p in points:
        d = dist(p, center)
        if d < radius:
            continue
        radius = .5 * (radius + d)
        old_to_new = d - radius
        new_center_x = (center[0] * radius + old_to_new * p[0]) / d
        new_center_y = (center[1] * radius + old_to_new * p[1]) / d
        new_center_z = (center[2] * radius + old_to_new * p[2]) / d
        center = (new_center_x, new_center_y, new_center_z)
    return center, radius
    
def generate_cube_points():
    x = uniform.rvs(size=(N, 3)).tolist()
    x[-1] = [0, 0, 0]
    x[-2] = [0, 0, 1]
    x[-3] = [1, 0, 0]
    x[-4] = [1, 1, 1]
    return x
    
def generate_sphere_points():
    def inside(c):
        return (c[0]-0.5)**2 + (c[1]-0.5)**2 + (c[2]-0.5)**2 < 1/4

    cube = uniform.rvs(size=(2*N, 3)).tolist()
    x = [c for c in cube if inside(c)]
    x[-1] = [0.5, 0.5, 0.0]
    x[-2] = [0.0, 0.5, 0.5]
    x[-3] = [0.5, 0.0, 0.5]
    x[-4] = [0.5, 1.0, 0.5]
    return x
   
def main():
    np.random.seed(seed=20210808)
    x = generate_cube_points()
    print("Cube")
    print(f"gist {bounding_sphere(x)}")
    print(f"ritter {bounding_sphere_ritter(x)}")
    

    x = generate_sphere_points()
    print("Sphere")
    print(f"gist {bounding_sphere(x)}")
    print(f"ritter {bounding_sphere_ritter(x)}")

    
if __name__ == "__main__":
    main()

ritter2.py

# Same code as above, with more runs to compare to the "exact" value.

import math
from scipy.stats import uniform
import numpy as np

def dist(a, b) -> float:
    return  ((a[0] - b[0])**2 + (a[1] - b[1])**2 + (a[2] - b[2])**2)**0.5
    
N = 10000
TRIALS = 10


def bounding_sphere(points):
    """
    From https://gist.github.com/paultsw/074b612d4e5980a18a94
    """
    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]):
            # Note that this step is flawed: bounding_sphere[0] (the center) never really changes)
            bounding_sphere = (bounding_sphere[0],dist(pt,bounding_sphere[0])) #grow the sphere to capture new point
    return bounding_sphere
   
def bounding_sphere_ritter(points):
    # Initial guess, same as before
    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
    center, radius = (((y[0]+z[0])/2,(y[1]+z[1])/2,(y[2]+z[2])/2), dist(y,z)/2)    #initial bounding sphere
    
    # Note this doesn't use the radius^2 optimization that Ritter uses
    for p in points:
        d = dist(p, center)
        if d < radius:
            continue
        radius = .5 * (radius + d)
        old_to_new = d - radius
        new_center_x = (center[0] * radius + old_to_new * p[0]) / d
        new_center_y = (center[1] * radius + old_to_new * p[1]) / d
        new_center_z = (center[2] * radius + old_to_new * p[2]) / d
        center = (new_center_x, new_center_y, new_center_z)
    return center, radius
    
def generate_cube_points():
    x = uniform.rvs(size=(N, 3)).tolist()
    x[-1] = [0, 0, 0]
    x[-2] = [0, 0, 1]
    x[-3] = [1, 0, 0]
    x[-4] = [1, 1, 1]
    return x
    
def generate_sphere_points():
    def inside(c):
        return (c[0]-0.5)**2 + (c[1]-0.5)**2 + (c[2]-0.5)**2 < 1/4

    cube = uniform.rvs(size=(2*N, 3)).tolist()
    x = [c for c in cube if inside(c)]
    x[-1] = [0.5, 0.5, 0.0]
    x[-2] = [0.0, 0.5, 0.5]
    x[-3] = [0.5, 0.0, 0.5]
    x[-4] = [0.5, 1.0, 0.5]
    return x
    
def run_trial(generator, exact):
    points = generator()
    c_gist, rad_gist = bounding_sphere(points)
    c_ritter, rad_ritter = bounding_sphere_ritter(points)
    
    # Make sure they're actually bounding spheres
    # Small safety margin to account for numerical precision
    for p in points:
        assert dist(p, c_gist) <= 1.000001 * rad_gist, f"{dist(p, c_gist)=}, {rad_gist=}"
        assert dist(p, c_ritter) <= 1.000001 * rad_ritter, f"{dist(p, c_ritter)=} {rad_ritter=}"
    error_gist = abs(rad_gist - exact)
    error_ritter = abs(rad_ritter - exact)
    return error_gist, error_ritter

def avg(x):
    return sum(x) / len(x)

def run_experiment(generator, exact):
    errors_gist = []
    errors_ritter = []
    for _ in range(TRIALS):
        error_gist, error_ritter = run_trial(generator, exact)
        errors_gist.append(error_gist)
        errors_ritter.append(error_ritter)
    print(f"gist  : max={max(errors_gist)}  avg={avg(errors_gist)}  pct={100.0 * max(errors_gist) / exact}%")
    print(f"ritter: max={max(errors_ritter)}  avg={avg(errors_ritter)}  pct={100.0 * max(errors_ritter) / exact}%")
   
def main():
    np.random.seed(seed=20210808)
    print("Cube")
    run_experiment(generate_cube_points, .5 * math.sqrt(3.0))
    
    print("Sphere")
    run_experiment(generate_sphere_points, .5)

    
if __name__ == "__main__":
    main()

Cube
gist ((0.5009538520928651, 0.48297642976169475, 0.49678676913174713), 0.8772593605335762)
ritter ((0.5072785458502712, 0.4895520489080106, 0.5031696262283336), 0.8661249358081996)

Sphere
gist ((0.48554103705593593, 0.5037648726351478, 0.5031219836853071), 0.5141840236538373)
ritter ((0.4917232591261705, 0.5013044493310901, 0.5024224321386422), 0.5076115940175479)

ritter2.py is the same calculations as the original, but run more times and compared against an "exact" answer.

The "bounding_sphere_ritter" approach is taken from Ritter's paper here https://www.researchgate.net/publication/242453691_An_Efficient_Bounding_Sphere

Output (values given are overestimates against the "exact" values):

gist  : max=0.018982604292868133  avg=0.010250889048852718  pct=2.1919223396815126%
ritter: max=0.005169140649164405  avg=0.0017581765408573925  pct=0.5968809490548213%
Sphere
gist  : max=0.04484039060824818  avg=0.027410374375479274  pct=8.968078121649636%
ritter: max=0.024441420887153065  avg=0.015042265461527437  pct=4.888284177430613%