improve volume estimation performance

core/post_processing.py

 """
 Post-processing utilities
+
+Notes: I don't have statistics on ranges of radii for each of the knots in
+the radial spline representation, but we regularly extract the average of
+these radii for each cell. So, depending on camera/lens, we get:
+    * 60x evolve: mean radii of 2-14 pixels (and measured areas of 30-750
+    pixels^2)
+    * 60x prime95b: mean radii of 3-24 pixels (and measured areas of 60-2000
+	pixels^2)
+
+And I presume that for a 100x lens we would get an ~5/3 increase over those
+values.
+
+In terms of the current volume estimation method, it's currently only
+implemented in the AnalysisToolbox repository, but it's super simple:
+
+mVol = 4/3*pi*sqrt(mArea/pi).^3
+
+where mArea is simply the sum of pixels for that cell.
 """
 import matplotlib.pyplot as plt
 import numpy as np
 from mpl_toolkits.mplot3d.art3d import Poly3DCollection
 from scipy import ndimage
 from skimage.morphology import erosion, ball
-from skimage import measure, draw, feature
+from skimage import measure, draw
+
+
+def my_ball(radius):
+    """Generates a ball-shaped structuring element.
+
+    This is the 3D equivalent of a disk.
+    A pixel is within the neighborhood if the Euclidean distance between
+    it and the origin is no greater than radius.
+
+    Parameters
+    ----------
+    radius : int
+        The radius of the ball-shaped structuring element.
+
+    Other Parameters
+    ----------------
+    dtype : data-type
+        The data type of the structuring element.
+
+    Returns
+    -------
+    selem : ndarray
+        The structuring element where elements of the neighborhood
+        are 1 and 0 otherwise.
+    """
+    n = 2 * radius + 1
+    Z, Y, X = np.mgrid[-radius:radius:n * 1j,
+              -radius:radius:n * 1j,
+              -radius:radius:n * 1j]
+    X **= 2
+    Y **= 2
+    Z **= 2
+    X += Y
+    X += Z
+    # s = X ** 2 + Y ** 2 + Z ** 2
+    return X <= radius * radius
+
+def circle_outline(r):
+    return ellipse_perimeter(r, r)
 
 def ellipse_perimeter(x, y):
     im_shape = int(2*max(x, y) + 1)
@@ -14,7 +71,7 @@ def ellipse_perimeter(x, y):
     rr, cc = draw.ellipse_perimeter(int(im_shape//2), int(im_shape//2),
                                     int(x), int(y))
     img[rr, cc] = 1
-    return img
+    return np.pad(img, 1)
 
 def capped_cylinder(x, y):
     max_size = (y + 2*x + 2)
@@ -34,48 +91,91 @@ def capped_cylinder(x, y):
     pixels ^= erosion(pixels)
     return pixels
 
+def volume_of_sphere(radius):
+    return 4 / 3 * np.pi * radius**3
+
+def plot_voxels(voxels):
+    verts, faces, normals, values = measure.marching_cubes_lewiner(
+        voxels, 0)
+    fig = plt.figure(figsize=(10, 10))
+    ax = fig.add_subplot(111, projection='3d')
+    mesh = Poly3DCollection(verts[faces])
+    mesh.set_edgecolor('k')
+    ax.add_collection3d(mesh)
+    ax.set_xlim(0, voxels.shape[0])
+    ax.set_ylim(0, voxels.shape[1])
+    ax.set_zlim(0, voxels.shape[2])
+    plt.tight_layout()
+    plt.show()
+
 # Volume estimation
-def union_of_spheres(outline, debug=False):
+def union_of_spheres(outline, shape='my_ball', debug=False):
     filled = ndimage.binary_fill_holes(outline)
     nearest_neighbor = ndimage.morphology.distance_transform_edt(
         outline == 0) * filled
     voxels = np.zeros((filled.shape[0], filled.shape[1], max(filled.shape)))
+    c_z = voxels.shape[2] // 2
     for x,y in zip(*np.where(filled)):
         radius = nearest_neighbor[(x,y)]
         if radius > 0:
-            b = ball(radius)
+            if shape == 'ball':
+                b = ball(radius)
+            elif shape == 'my_ball':
+                b = my_ball(radius)
+            else:
+                raise ValueError(f"{shape} is not an accepted value for "
+                                 f"shape.")
             centre_b = ndimage.measurements.center_of_mass(b)
 
             I,J,K = np.ogrid[:b.shape[0], :b.shape[1], :b.shape[2]]
-            c_z = voxels.shape[2]//2
             voxels[I + int(x - centre_b[0]), J + int(y - centre_b[1]),
                    K + int(c_z - centre_b[2])] += b
     if debug:
-        verts, faces, normals, values = measure.marching_cubes_lewiner(
-            voxels, 0)
-        fig = plt.figure(figsize=(10, 10))
-        ax = fig.add_subplot(111, projection='3d')
-        mesh = Poly3DCollection(verts[faces])
-        mesh.set_edgecolor('k')
-        ax.add_collection3d(mesh)
-        ax.set_xlim(0, filled.shape[0])
-        ax.set_ylim(0, filled.shape[1])
-        ax.set_zlim(0, max(filled.shape))
-        plt.tight_layout()
-        plt.show()
+        plot_voxels(voxels)
     return voxels.astype(bool).sum()
 
-def conical(outline, debug=False, selem=None):
+def improved_uos(outline, shape='my_ball', debug=False):
     filled = ndimage.binary_fill_holes(outline)
-    cone = [filled]
-    while filled.sum() > 0 :
-        filled = erosion(filled, selem=selem)
-        cone.append(filled)
-        if debug:
-            plt.imshow(filled)
-            plt.show()
-    cone = np.dstack(cone)
-    return 4 * np.sum(cone) #* 0.95 #To make the circular version work
+    nearest_neighbor = ndimage.morphology.distance_transform_edt(
+        outline == 0) * filled
+    voxels = np.zeros((filled.shape[0], filled.shape[1], max(filled.shape)))
+    c_z = voxels.shape[2] // 2
+
+    while np.any(nearest_neighbor != 0):
+        radius = np.max(nearest_neighbor)
+        x, y = np.argwhere(nearest_neighbor == radius)[0]
+        if shape == 'ball':
+            b = ball(np.ceil(radius))
+        elif shape == 'my_ball':
+            b = my_ball(np.ceil(radius))
+        else:
+            raise ValueError(f"{shape} is not an accepted value for shape")
+        centre_b = ndimage.measurements.center_of_mass(b)
+
+        I, J, K = np.ogrid[:b.shape[0], :b.shape[1], :b.shape[2]]
+        voxels[I + int(x - centre_b[0]), J + int(y - centre_b[1]),
+               K + int(c_z - centre_b[2])] += b
+
+        # Use the central disk of the ball from voxels to get the circle
+        # = 0 if nn[x,y] < r else nn[x,y]
+        rr, cc = draw.circle(x, y, np.ceil(radius), nearest_neighbor.shape)
+        nearest_neighbor[rr, cc] = 0
+    if debug:
+        plot_voxels(voxels)
+    return voxels.astype(bool).sum()
+
+def conical(outline, debug=False):
+    nearest_neighbor = ndimage.morphology.distance_transform_edt(
+        outline == 0) * ndimage.binary_fill_holes(outline)
+    if debug:
+        hf = plt.figure()
+        ha = hf.add_subplot(111, projection='3d')
+
+        X, Y = np.meshgrid(np.arange(nearest_neighbor.shape[0]),
+                           np.arange(nearest_neighbor.shape[1]))
+        ha.plot_surface(X, Y, nearest_neighbor)
+        plt.show()
+    return 4 * nearest_neighbor.sum()
 
 def volume(outline, method='spheres'):
     if method=='conical':

test/test_post_processing.py

 import matplotlib.pyplot as plt
 import numpy as np
 import skimage.morphology as morph
+from scipy import ndimage
 from skimage import draw
 import unittest
 
-from core.post_processing import conical, ellipse_perimeter, union_of_spheres
+from core.post_processing import conical, ellipse_perimeter, \
+    union_of_spheres, volume_of_sphere, circle_outline
 
 
 class VolumeEstimation(unittest.TestCase):
     def test_conical(self):
         radius = np.random.choice(range(60, 100))
-        con = conical(morph.disk(radius))
+        con = conical(circle_outline(radius))
         b_sum = morph.ball(radius).sum()
         # Close relative to the value.
         print(radius, con, b_sum)
@@ -19,13 +21,13 @@ class VolumeEstimation(unittest.TestCase):
     def test_conical_ellipse(self):
         e = ellipse_perimeter(4, 5)
         con = conical(e)
-        true = draw.ellipsoid_stats(4, 4, 5)[0]
+        true = draw.ellipsoid_stats(4, 5, 4)[0]
         print(con, true)
 
     def test_sphere_error(self):
         radii = range(3, 30)
-        con = [conical(morph.disk(radius)) for radius in radii]
-        spheres = [union_of_spheres(ellipse_perimeter(r, r)) for r in radii]
+        con = [conical(circle_outline(radius)) for radius in radii]
+        spheres = [union_of_spheres(circle_outline(r)) for r in radii]
         true = [4*(r**3)*np.pi/3 for r in radii]
         mVol = [4 / 3 * np.pi * np.sqrt(morph.disk(radius).sum() / np.pi)**3
                 for radius in radii]
@@ -42,15 +44,15 @@ class VolumeEstimation(unittest.TestCase):
 
     def test_ellipse_error(self):
         x_radii = range(3, 30)
-        y_radii = [np.ceil(1.2*r) for r in x_radii]
+        y_radii = [np.ceil(2.5*r) for r in x_radii]
         ellipses = [ellipse_perimeter(x_r, y_r)
                     for x_r, y_r in zip(x_radii, y_radii)]
         con = [conical(ellipse) for ellipse in ellipses]
         spheres = [union_of_spheres(ellipse) for ellipse in ellipses]
-        mVol = [(4 * np.pi * np.sqrt(ellipse.sum() / np.pi)**3) / 3
-                for ellipse in ellipses]
-        true = [draw.ellipsoid_stats(x_r, y_r, x_r)[0]
-                for x_r, y_r in zip(x_radii, y_radii)]
+        mVol = np.array([4 / 3 * np.pi * np.sqrt(ndimage.binary_fill_holes(
+            ellipse).sum() / np.pi) ** 3 for ellipse in ellipses])
+        true = np.array([4 * np.pi * x_r * y_r * x_r / 3
+                         for x_r, y_r in zip(x_radii, y_radii)])
         plt.scatter(true, con, label='Conical')
         plt.scatter(true, spheres, label='Spheres')
         plt.scatter(true, mVol, label='mVol')
@@ -61,8 +63,32 @@ class VolumeEstimation(unittest.TestCase):
         plt.legend()
         plt.show()
 
-    def test_mixed_error(self):
-        pass
+    def test_minor_major_error(self):
+        r = np.random.choice(list(range(3, 30)))
+        x_radii = np.linspace(r/3, r, 20)
+        y_radii = r**2/x_radii
+
+        ellipses = [ellipse_perimeter(x_r, y_r)
+                    for x_r, y_r in zip(x_radii, y_radii)]
+        con = np.array([conical(ellipse) for ellipse in ellipses])
+        spheres = np.array([union_of_spheres(ellipse) for ellipse in ellipses])
+        mVol = np.array([4 / 3 * np.pi * np.sqrt(ndimage.binary_fill_holes(
+            ellipse).sum() / np.pi) ** 3
+                for ellipse in ellipses])
+
+        true = np.array([4 * np.pi * x_r * y_r * x_r / 3
+                         for x_r, y_r in zip(x_radii, y_radii)])
+
+        ratio = y_radii/x_radii
+        plt.scatter(ratio, con/true, label='Conical')
+        plt.scatter(ratio, spheres/true, label='Spheres')
+        plt.scatter(ratio, mVol/true, label='mVol')
+        plt.xlabel("Major/Minor")
+        plt.ylabel("Estimated / Analytical")
+        plt.title(f"Error by circularity, r = {r}")
+        plt.legend()
+        plt.show()
+