docs: 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 []))