compute ellipticities of a light distribution

    :param I_xy:
    :param x:
    :param y:
    :return:
192172
def ellipticities(I_xy, x, y):
    
    Q_xx, Q_xy, Q_yy, bkg = moments(I_xy, x, y)
    norm = Q_xx + Q_yy + 2 * np.sqrt(Q_xx*Q_yy - Q_xy**2)
    e1 = (Q_xx - Q_yy) / norm
    e2 = 2 * Q_xy / norm
    return e1 / (1+bkg), e2 / (1+bkg)
897043
def _2ndderiv_xyz(self,x,y,z,i,j):
        """"""
        return -np.pi*self._rhoc_M*self.a**3*self._b*self._c *\
            _2ndDerivInt(x,y,z,self._a2,self._b2*self._a2,self._c2*self._a2,self.n,i,j)
386902
def pixel(self, ix, iy):
        
        x = (ix - self.xoffset + 0.5) / self.scale
        y = (iy - self.yoffset + 0.5) / self.scale
        z = complex(x, y)
        self.set_c(z)
        return self.mpixel(z)
1137745
def get_gradient_y(shape, py):
    
    import scipy.sparse

    height, width = shape
    rows = []
    empty = scipy.sparse.dia_matrix((width, width))
    identity = scipy.sparse.identity(width)

    # Create the blocks by row, beginning with blocks leading up to the peak row from the top
    for n in range(py):
        row = [empty]*n
        row += [-identity, identity]
        row += [empty]*(height-n-2)
        rows.append(row)
    # Set all elements in the peak row to zero
    rows.append([empty]*height)
    # Create the blocks for the rows leading up to the peak row from the bottom
    for n in range(height-py-1):
        row = [empty]*(py+n)
        row += [identity, -identity]
        row += [empty]*(height-py-n-2)
        rows.append(row)
    return scipy.sparse.bmat(rows)
470981
def in_lamp_reach(p):
    '''  '''
    v1 = XYPoint(Lime.x - Red.x, Lime.y - Red.y)
    v2 = XYPoint(Blue.x - Red.x, Blue.y - Red.y)

    q = XYPoint(p.x - Red.x, p.y - Red.y)
    s = cross_product(q, v2) / cross_product(v1, v2)
    t = cross_product(v1, q) / cross_product(v1, v2)

    return (s >= 0.0) and (t >= 0.0) and (s + t <= 1.0)
405274
def get_xy_1D(ds, stride=1, getval=False):
    
    gt = ds.GetGeoTransform()
    #stride = stride_m/gt[1]
    pX = np.arange(0, ds.RasterXSize, stride)
    pY = np.arange(0, ds.RasterYSize, stride)
    mX, dummy = pixelToMap(pX, pY[0], gt)
    dummy, mY = pixelToMap(pX[0], pY, gt)
    return mX, mY
77985
def shakeshake2_py(x, y, equal=False, individual=False):
  """"""
  if equal:
    alpha = 0.5
  elif individual:
    alpha = tf.random_uniform(tf.get_shape(x)[:1])
  else:
    alpha = tf.random_uniform([])

  return alpha * x + (1.0 - alpha) * y
244980
def DSYR_blas(A, x, alpha=1.):
    
    At = blas.dsyr(lower=0, x=x, a=A, alpha=alpha, overwrite_a=False) #See https://github.com/scipy/scipy/issues/8155
    A[:] = At
    symmetrify(A, upper=True)
1349410
def _xys(date):
    

    X, Y, s_xy2 = _xysxy2(date)

    # convert milli-arcsecond to arcsecond
    dX, dY = date.eop.dx / 1000., date.eop.dy / 1000.

    # Convert arcsecond to degrees then to radians
    X = np.radians((X + dX) / 3600.)
    Y = np.radians((Y + dY) / 3600.)
    s = np.radians(s_xy2 / 3600.) - (X * Y / 2)

    return X, Y, s
1358952
static public double xy(double[] x, double[] y) {
    double result = 0.0;
    for (int i = 0; i < x.length; i++) {
      result += (x[i] * y[i]);
    }

    return (result);
  }
1335550
def example_ylm(m=0, n=2, shape=128, limits=[-4, 4], draw=True, show=True, **kwargs):
    """"""
    import ipyvolume.pylab as p3

    __, __, __, r, theta, phi = xyz(shape=shape, limits=limits, spherical=True)
    radial = np.exp(-(r - 2) ** 2)
    data = np.abs(scipy.special.sph_harm(m, n, theta, phi) ** 2) * radial  # pylint: disable=no-member
    if draw:
        vol = p3.volshow(data=data, **kwargs)
        if show:
            p3.show()
        return vol
    else:
        return data
1222606
def xyz2rgb(x__, y__, z__):
    
    x2_ = x__ / 100.0
    y2_ = y__ / 100.0
    z2_ = z__ / 100.0

    r__ = x2_ *  3.2406 + y2_ * -1.5372 + z2_ * -0.4986
    g__ = x2_ * -0.9689 + y2_ *  1.8758 + z2_ *  0.0415
    b__ = x2_ *  0.0557 + y2_ * -0.2040 + z2_ *  1.0570

    def finv(arr):
        """Inverse
        """
        return np.where(arr > 0.0031308,
                        1.055 * (arr ** (1.0 / 2.4)) - 0.055,
                        12.92 * arr)

    return finv(r__) * 255, finv(g__) * 255, finv(b__) * 255
443947
def C(w, Xs):
    '''
    '''
    n = len(Xs)
    P = projection_matrix(w)
    Ys = [np.dot(P, X) for X in Xs]
    A = calc_A(Ys)
    A_hat = calc_A_hat(A, skew_matrix(w))

    return np.dot(A_hat, sum(np.dot(Y, Y) * Y for Y in Ys)) / np.trace(np.dot(A_hat, A))
1031681
def inROI(self, Y):
        ''''''
        if Y.ndim > 1:
            area = np.zeros((Y.shape[0],4))
        else:
            area = np.zeros((1,4))
        pts = np.zeros((0,), int)
        pdist = np.zeros((0,), int)
        dist0 = 0
        for k in range(len(self.prect)):
            self.square_area = (triangle_area(self.prect[k][0,:], self.prect[k][1,:], self.prect[k][2,:]) +
                            triangle_area(self.prect[k][2,:], self.prect[k][3,:], self.prect[k][4,:]))
            for n in range(4):
                area[:,n] = triangle_area(self.prect[k][0+n,:], self.prect[k][1+n,:], Y)
            # points inside prect
            newpts = np.array((area.sum(axis=1) <= self.square_area+1e-5).nonzero()).flatten().astype(int)
            if newpts.size > 0:
                pts = np.concatenate((pts, newpts))
                newdists = self.orthproj(Y[newpts, :], k) + dist0
                pdist = np.concatenate((pdist, newdists))
            dist0 += (np.diff(self.pos[k], axis=0)[0,:]**2).sum()
            # check if in radius of circle
            if k < len(self.prect)-1:
                pcent = self.pos[k][1,:]
                dist = ((Y - pcent[np.newaxis,:])**2).sum(axis=1)**0.5
                newpts = np.array((dist<=self.d).nonzero()[0].astype(int))
                if newpts.size > 0:
                    pts = np.concatenate((pts, newpts))
                    newdists = dist0 * np.ones(newpts.shape)
                    pdist = np.concatenate((pdist, newdists))

        pts, inds = np.unique(pts, return_index=True)
        pdist = pdist[inds]
        return pts, pdist
375536
private void computePsiSmooth(GrayF32 ux , GrayF32 uy , GrayF32 vx , GrayF32 vy ,
								  GrayF32 psiSmooth ) {
		int N = derivFlowUX.width * derivFlowUX.height;

		for( int i = 0; i < N; i++ ) {
			float vux = ux.data[i];
			float vuy = uy.data[i];
			float vvx = vx.data[i];
			float vvy = vy.data[i];

			float mu = vux*vux + vuy*vuy;
			float mv = vvx*vvx + vvy*vvy;

			psiSmooth.data[i] = (float)(1.0/(2.0*Math.sqrt(mu + mv + EPSILON*EPSILON)));
		}
	}
1219957
protected float intersectionY (IRay3 ray, float y) {
        IVector3 origin = ray.origin(), dir = ray.direction();
        float t = (y - origin.y()) / dir.y();
        if (t < 0f) {
            return Float.MAX_VALUE;
        }
        float ix = origin.x() + t*dir.x(), iz = origin.z() + t*dir.z();
        return (ix >= _minExtent.x && ix <= _maxExtent.x &&
            iz >= _minExtent.z && iz <= _maxExtent.z) ? t : Float.MAX_VALUE;
    }
743362
def gaussian2d(xy, sx, sy, mx=0, my=0, rho=0, amp=1, offs=0):
    '''
    
    '''
    x,y = xy
    return offs+amp*( 
        1/(2*np.pi*sx*sy*(1-(rho**2))**0.5) *
         np.exp( (-1/(2*(1-rho**2))) *
                 (
                    ( (x-mx)**2/sx**2 )
                  + ( (y-my)**2/sy**2 )
                  - ( ( 2*rho*(x-mx)*(y-my)) / (sx*sy) )
                  )
                )
         )