--- a
+++ b/CLI/MusculoskeletalAnalysisCLITools/thickness.py
@@ -0,0 +1,41 @@
+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))