def findSpheres(mask):
    """Calculates thickness by finding the largest sphere containing each point."""
    import numpy as np
    from math import floor
    from scipy.ndimage import distance_transform_edt

    dist = distance_transform_edt(mask)
    rads = dist.copy()

    a,b,c=np.nonzero(mask)
    order = np.argsort(dist[a,b,c])
    for n in order:
        # For each point, find the thickness radius at that point from the distance map
        x=a[n]
        y=b[n]
        z=c[n]

        r=dist[x,y,z]
        rf = floor(r)
        # Check that ranges are inside boundaries
        xr0=min(rf,x)
        xr1=min(rf, dist.shape[0]-x-1)
        yr0=min(rf,y)
        yr1=min(rf, dist.shape[1]-y-1)
        zr0=min(rf,z)
        zr1=min(rf, dist.shape[2]-z-1)
        # Creates a cubic view of the region of interest
        roi=rads[x-xr0:x+xr1+1, y-yr0:y+yr1+1, z-zr0:z+zr1+1]
        # Creates mask of sphere within radius
        i,j,k=np.mgrid[-xr0:xr1+1,-yr0:yr1+1,-zr0:zr1+1]
        mask = approxLTE(i**2 + j**2 + k**2, np.full_like(roi, r**2)) 
        mask = np.bitwise_and(mask, roi>0)
        # Sets all values inside region of interest and sphere mask to the radius if it is higher than their current value
        mask = mask * r
        roi[:,:,:] = np.maximum(roi, mask)
    return rads[np.nonzero(rads)]

def approxLTE(a, b):
    # Peforms array elementwise less than or equal comparisons with some tolerance for floating point values
    import numpy as np
    return np.bitwise_or(a <= b, np.isclose(a,b))