code: def moments(I_xy_input, x, y):
      
      I_xy = copy.deepcopy(I_xy_input)
      background = np.minimum(0, np.min(I_xy))
      I_xy -= background
      x_ = np.sum(I_xy * x)
      y_ = np.sum(I_xy * y)
      r = (np.max(x) - np.min(x)) / 3.
      mask = mask_util.mask_sphere(x, y, center_x=x_, center_y=y_, r=r)
      Q_xx = np.sum(I_xy * mask * (x - x_) ** 2)
      Q_xy = np.sum(I_xy * mask * (x - x_) * (y - y_))
      Q_yy = np.sum(I_xy * mask * (y - y_) ** 2)
      return Q_xx, Q_xy, Q_yy, background / np.mean(I_xy)
docs: compute quadrupole moments from a light distribution
  
      :param I_xy: light distribution
      :param x: x-coordinates of I_xy
      :param y: y-coordinates of I_xy
      :return: Q_xx, Q_xy, Q_yy