ellipsoid.py

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