def ellipse_from_second_moments_ijv(i,j, image, labels, indexes, wants_compactness = False):
    """Calculate measurements of ellipses equivalent to the second moments of labels
    
    i,j - coordinates of each point
    image  - the intensity at each point
    labels - for each labeled object, derive an ellipse
    indexes - sequence of indexes to process
    
    returns the following arrays:
       coordinates of the center of the ellipse
       eccentricity
       major axis length
       minor axis length
       orientation
    
    some definitions taken from "Image Moments-Based Structuring and Tracking
    of Objects", LOURENA ROCHA, LUIZ VELHO, PAULO CEZAR P. CARVALHO,
    http://sibgrapi.sid.inpe.br/col/sid.inpe.br/banon/2002/10.23.11.34/doc/35.pdf
    particularly equation 5 (which has some errors in it).
    These yield the rectangle with equivalent second moments. I translate
    to the ellipse by multiplying by 1.154701 which is Matlab's calculation
    of the major and minor axis length for a square of length X divided
    by the actual length of the side of a square of that length.
    
    eccentricity is the distance between foci divided by the major axis length
    orientation is the angle of the major axis with respect to the X axis
    """
    if len(indexes) == 0:
        return [np.zeros((0,2))] + [np.zeros((0,))] * (5 if wants_compactness else 4)
    if len(i) == 0:
        return ([np.zeros((len(indexes), 2)), np.ones(len(indexes))] +
                [np.zeros(len(indexes))] * (4 if wants_compactness else 3))
    #
    # Normalize to center of object for stability
    #
    nlabels = np.max(indexes)+1
    m = np.array([[None, 0, None],
                  [0, None, None],
                  [None, None, None]], object)
    if np.all(image == 1):
        image = 1
        m[0,0] = intensity = np.bincount(labels)
    else:
        m[0,0] = intensity = np.bincount(labels, image)
    ic = np.bincount(labels, i * image) / intensity
    jc = np.bincount(labels, j * image) / intensity
    i = i - ic[labels]
    j = j - jc[labels]
    #
    # Start by calculating the moments m[p][q] of the image
    # sum(i**p j**q)
    #
    # m[1,0] = 0 via normalization
    # m[0,1] = 0 via normalization
    m[1,1] = np.bincount(labels, i*j*image)
    m[2,0] = np.bincount(labels, i*i*image)
    m[0,2] = np.bincount(labels, j*j*image)
    
    a = m[2,0] / m[0,0]
    b = 2*m[1,1]/m[0,0]
    c = m[0,2] / m[0,0]
    
    theta = np.arctan2(b,c-a) / 2
    temp = np.sqrt(b**2+(a-c)**2)
    #
    # If you do a linear regression of the circles from 1 to 50 radius
    # in Matlab, the resultant values fit a line with slope=.9975 and
    # intercept .095. I'm adjusting the lengths accordingly.
    #
    mystery_constant = 0.095
    mystery_multiplier = 0.9975
    major_axis_len = (np.sqrt(8*(a+c+temp)) * mystery_multiplier +
                      mystery_constant)
    minor_axis_len = (np.sqrt(8*(a+c-temp)) * mystery_multiplier +
                      mystery_constant)
    eccentricity = np.sqrt(1-(minor_axis_len / major_axis_len)**2)
    compactness = 2 * np.pi * (a + c) / m[0,0]
    return ([np.column_stack((ic[indexes], jc[indexes])),
             eccentricity[indexes],
             major_axis_len[indexes],
             minor_axis_len[indexes],
             theta[indexes]] +
            ([compactness[indexes]] if wants_compactness else []))