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