def angular_distribution(labels, resolution=100, weights=None):
    '''For each object in labels, compute the angular distribution
    around the centers of mass.  Returns an i x j matrix, where i is
    the number of objects in the label matrix, and j is the resolution
    of the distribution (default 100), mapped from -pi to pi.

    Optionally, the distributions can be weighted by pixel.

    The algorithm approximates the angular width of pixels relative to
    the object centers, in an attempt to be accurate for small
    objects.

    The ChordRatio of an object can be approximated by 
    >>> angdist = angular_distribution(labels, resolution)
    >>> angdist2 = angdist[:, :resolution//2] + angdist[:, resolution//2] # map to widths, rather than radii
    >>> chord_ratio = np.sqrt(angdist2.max(axis=1) / angdist2.min(axis=1)) # sqrt because these are sectors, not triangles
    '''
    if weights is None:
        weights = np.ones(labels.shape)
    maxlabel = labels.max()
    ci, cj = centers_of_labels(labels)
    j, i = np.meshgrid(np.arange(labels.shape[0]), np.arange(labels.shape[1]))
    # compute deltas from pixels to object centroids, and mask to labels
    di = i[labels > 0] - ci[labels[labels > 0] - 1]
    dj = j[labels > 0] - cj[labels[labels > 0] - 1]
    weights = weights[labels > 0]
    labels = labels[labels > 0]
    # find angles, and angular width of pixels
    angle = np.arctan2(di, dj)
    # Use pixels of width 2 to get some smoothing
    width = np.arctan(1.0 / np.sqrt(di**2 + dj**2 + np.finfo(float).eps))
    # create an onset/offset array of size 3 * resolution
    lo = np.clip((angle - width) * resolution / (2 * np.pi), -resolution, 2 * resolution).astype(int) + resolution
    hi = np.clip((angle + width) * resolution / (2 * np.pi), -resolution, 2 * resolution).astype(int) + resolution
    # make sure every pixel counts at least once
    hi[lo == hi] += 1
    # normalize weights by their angular width (adding a bit to avoid 0 / 0)
    weights /= (hi - lo)
    onset = scipy.sparse.coo_matrix((weights, (labels - 1, lo)), (maxlabel, 3 * resolution)).toarray()
    offset = scipy.sparse.coo_matrix((weights, (labels - 1, hi)), (maxlabel, 3 * resolution)).toarray()
    # sum onset/offset to get actual distribution
    onoff = onset - offset
    dist = np.cumsum(onoff, axis=1)
    dist = dist[:, :resolution] + dist[:, resolution:2*resolution] + dist[:, 2*resolution:]
    return dist