dgobbi · November 21, 2017 20:45 · Nov 21, 2017
diff --git a/ellipsoid.py b/ellipsoid.py
@@ -0,0 +1,35 @@
+import numpy as np
+
+def ellipsoid(radii, dtype='d'):
+    """Generate an ellipsoid with the specified radii.
+
+    The ellipsoid will be centered in the output, which will always
+    have odd dimensions.  The surface of the sphere will be antialiased,
+    that is, a linear ramp is applied at the surface instead of a hard
+    edge.  The size of the ramp is exactly one sample spacing, i.e. points
+    right on the surface will have a value of 0.5 while points that are
+    half a sample spacing or more from the edge, measured along the surface
+    normal, will have a value of 1.0 (inside) or 0.0 (outside).
+    """
+    # the dimensions will always be odd.
+    shape = [2*int(np.ceil(r))+1 for r in radii]
+    grid = []
+    flat_center = 0
+    for m,r in zip(shape, radii):
+        x = np.mgrid[0:m].astype('d')
+        x -= 0.5*(m-1)
+        x *= 1.0/r
+        grid.append(x)
+        flat_center = m*flat_center + m//2
+    mg = np.meshgrid(*grid,indexing='ij')
+    # compute radius from voxel position
+    distance = np.sqrt(sum(np.square(a) for a in mg))
+    distance.flat[flat_center] = 1.0
+    # approximate signed distance from the ellipsoid
+    vec = [x/distance*r for x,r in zip(mg,radii)]
+    distance -= 1.0
+    distance *= np.sqrt(sum(np.square(a) for a in vec))
+    distance.flat[flat_center] = -min(radii)
+    # generate a mask, use ramp to avoid discontinuity
+    mask = np.interp(distance, [-0.5,+0.5], [1.0, 0.0])
+    return mask