Volume estimation first tests.

core/post_processing.py

+"""
+Post-processing utilities
+"""
+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
+
+def ellipse_perimeter(x, y):
+    im_shape = int(2*max(x, y) + 1)
+    img = np.zeros((im_shape, im_shape), dtype=np.uint8)
+    rr, cc = draw.ellipse_perimeter(int(im_shape//2), int(im_shape//2),
+                                    int(x), int(y))
+    img[rr, cc] = 1
+    return img
+
+def capped_cylinder(x, y):
+    max_size = (y + 2*x + 2)
+    pixels = np.zeros((max_size, max_size))
+
+    rect_start = ((max_size-x)//2, x + 1)
+    rr, cc = draw.rectangle_perimeter(rect_start, extent=(x, y),
+                                     shape=(max_size, max_size))
+    pixels[rr, cc] = 1
+    circle_centres = [(max_size//2 - 1, x),
+                      (max_size//2 - 1, max_size - x - 1 )]
+    for r, c in circle_centres:
+        rr, cc = draw.circle_perimeter(r, c, (x + 1)//2,
+                                       shape=(max_size, max_size))
+        pixels[rr, cc] = 1
+    pixels = ndimage.morphology.binary_fill_holes(pixels)
+    pixels ^= erosion(pixels)
+    return pixels
+
+# Volume estimation
+def union_of_spheres(outline, 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)))
+    for x,y in zip(*np.where(filled)):
+        radius = nearest_neighbor[(x,y)]
+        if radius > 0:
+            b = ball(radius)
+            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()
+    return voxels.astype(bool).sum()
+
+def conical(outline, debug=False, selem=None):
+    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
+
+def volume(outline, method='spheres'):
+    if method=='conical':
+        return conical(outline)
+    elif method=='spheres':
+        return union_of_spheres(outline)
+    else:
+        raise ValueError(f"Method {method} not implemented.")
+
+def circularity(outline):
+    pass
\ No newline at end of file

core/utils.py

@@ -44,5 +44,5 @@ class ImageCache:
         # Clean up the queue
         self._queue.append(item)
         if len(self._queue) > self.max_len:
-            del self._dict.pop[self._queue.pop(0)]
+            del self._dict[self._queue.pop(0)]
 

test/test_post_processing.py

+import matplotlib.pyplot as plt
+import numpy as np
+import skimage.morphology as morph
+from skimage import draw
+import unittest
+
+from core.post_processing import conical, ellipse_perimeter, union_of_spheres
+
+
+class VolumeEstimation(unittest.TestCase):
+    def test_conical(self):
+        radius = np.random.choice(range(60, 100))
+        con = conical(morph.disk(radius))
+        b_sum = morph.ball(radius).sum()
+        # Close relative to the value.
+        print(radius, con, b_sum)
+        self.assertAlmostEqual(abs(con - b_sum)/b_sum, 0, delta=0.10)
+
+    def test_conical_ellipse(self):
+        e = ellipse_perimeter(4, 5)
+        con = conical(e)
+        true = draw.ellipsoid_stats(4, 4, 5)[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]
+        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]
+        plt.scatter(true, con, label='Conical')
+        plt.scatter(true, spheres, label='Spheres')
+        plt.scatter(true, mVol, label='mVol')
+        plt.plot(true, true, 'k-' )
+        plt.xlabel("Analytical")
+        plt.ylabel("Estimated")
+        plt.title("Disk")
+        plt.legend()
+        plt.show()
+
+
+    def test_ellipse_error(self):
+        x_radii = range(3, 30)
+        y_radii = [np.ceil(1.2*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)]
+        plt.scatter(true, con, label='Conical')
+        plt.scatter(true, spheres, label='Spheres')
+        plt.scatter(true, mVol, label='mVol')
+        plt.plot(true, true, 'k-')
+        plt.xlabel("Analytical")
+        plt.ylabel("Estimated")
+        plt.title("Ellipse")
+        plt.legend()
+        plt.show()
+
+    def test_mixed_error(self):
+        pass
+
+
+
+
+
+if __name__ == '__main__':
+    unittest.main()