Python numpy.real() 使用实例

The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.

Example 1

def encode(img_path, wm_path, res_path, alpha):
    img = cv2.imread(img_path)
    img_f = np.fft.fft2(img)
    height, width, channel = np.shape(img)
    watermark = cv2.imread(wm_path)
    wm_height, wm_width = watermark.shape[0], watermark.shape[1]
    x, y = range(height / 2), range(width)
    random.seed(height + width)
    random.shuffle(x)
    random.shuffle(y)
    tmp = np.zeros(img.shape)
    for i in range(height / 2):
        for j in range(width):
            if x[i] < wm_height and y[j] < wm_width:
                tmp[i][j] = watermark[x[i]][y[j]]
                tmp[height - 1 - i][width - 1 - j] = tmp[i][j]
    res_f = img_f + alpha * tmp
    res = np.fft.ifft2(res_f)
    res = np.real(res)
    cv2.imwrite(res_path, res, [int(cv2.IMWRITE_JPEG_QUALITY), 100]) 

Example 2

def _ncc_c(x, y):
    """
    >>> _ncc_c([1,2,3,4], [1,2,3,4])
    array([ 0.13333333,  0.36666667,  0.66666667,  1.        ,  0.66666667,
            0.36666667,  0.13333333])
    >>> _ncc_c([1,1,1], [1,1,1])
    array([ 0.33333333,  0.66666667,  1.        ,  0.66666667,  0.33333333])
    >>> _ncc_c([1,2,3], [-1,-1,-1])
    array([-0.15430335, -0.46291005, -0.9258201 , -0.77151675, -0.46291005])
    """
    den = np.array(norm(x) * norm(y))
    den[den == 0] = np.Inf

    x_len = len(x)
    fft_size = 1<<(2*x_len-1).bit_length()
    cc = ifft(fft(x, fft_size) * np.conj(fft(y, fft_size)))
    cc = np.concatenate((cc[-(x_len-1):], cc[:x_len]))
    return np.real(cc) / den 

Example 3

def quaternion_real(quaternion):
    """Return real part of quaternion.

    >>> quaternion_real([3, 0, 1, 2])
    3.0

    """
    return float(quaternion[0]) 

Example 4

def correlate(self, imgfft):
        #Very much related to the convolution theorem, the cross-correlation
        #theorem states that the Fourier transform of the cross-correlation of
        #two functions is equal to the product of the individual Fourier
        #transforms, where one of them has been complex conjugated:  
        
        
        if self.imgfft is not 0 or imgfft.imgfft is not 0:
            imgcj = np.conjugate(self.imgfft)
            imgft = imgfft.imgfft
            
            prod = deepcopy(imgcj)
            for x in range(imgcj.shape[0]):
                for y in range(imgcj.shape[0]):
                    prod[x][y] = imgcj[x][y] * imgft[x][y]
            
            cc = Corr( np.real(fft.ifft2(fft.fftshift(prod)))) # real image of the correlation
            
            # adjust to center
            cc.data = np.roll(cc.data, int(cc.data.shape[0] / 2), axis = 0)
            cc.data = np.roll(cc.data, int(cc.data.shape[1] / 2), axis = 1)
        else:
            raise FFTnotInit()
        return cc 

Example 5

def save_image(imager, grid_data, grid_norm, output_file):
    """Makes an image from gridded visibilities and saves it to a FITS file.

    Args:
        imager (oskar.Imager):          Handle to configured imager.
        grid_data (numpy.ndarray):      Final visibility grid.
        grid_norm (float):              Grid normalisation to apply.
        output_file (str):              Name of output FITS file to write.
    """
    # Make the image (take the FFT, normalise, and apply grid correction).
    imager.finalise_plane(grid_data, grid_norm)
    grid_data = numpy.real(grid_data)

    # Trim the image if required.
    border = (imager.plane_size - imager.image_size) // 2
    if border > 0:
        end = border + imager.image_size
        grid_data = grid_data[border:end, border:end]

    # Write the FITS file.
    hdr = fits.header.Header()
    fits.writeto(output_file, grid_data, hdr, clobber=True) 

Example 6

def nufft_scale1(N, K, alpha, beta, Nmid):
    '''
    calculate image space scaling factor
    '''
#     import types
#     if alpha is types.ComplexType:
    alpha = numpy.real(alpha)
#         print('complex alpha may not work, but I just let it as')

    L = len(alpha) - 1
    if L > 0:
        sn = numpy.zeros((N, 1))
        n = numpy.arange(0, N).reshape((N, 1), order='F')
        i_gam_n_n0 = 1j * (2 * numpy.pi / K) * (n - Nmid) * beta
        for l1 in range(-L, L + 1):
            alf = alpha[abs(l1)]
            if l1 < 0:
                alf = numpy.conj(alf)
            sn = sn + alf * numpy.exp(i_gam_n_n0 * l1)
    else:
        sn = numpy.dot(alpha, numpy.ones((N, 1), dtype=numpy.float32))
    return sn 

Example 7

def kaiser_bessel_ft(u, J, alpha, kb_m, d):
    '''
    Interpolation weight for given J/alpha/kb-m
    '''

    u = u * (1.0 + 0.0j)
    import scipy.special
    z = numpy.sqrt((2 * numpy.pi * (J / 2) * u) ** 2.0 - alpha ** 2.0)
    nu = d / 2 + kb_m
    y = ((2 * numpy.pi) ** (d / 2)) * ((J / 2) ** d) * (alpha ** kb_m) / \
        scipy.special.iv(kb_m, alpha) * scipy.special.jv(nu, z) / (z ** nu)
    y = numpy.real(y)
    return y 

Example 8

def mdct(x, odd=True):
    """ Calculate modified discrete cosine transform of input signal

    Parameters
    ----------
    X : array_like
        The input signal
    odd : boolean, optional
        Switch to oddly stacked transform. Defaults to :code:`True`.

    Returns
    -------
    out : array_like
        The output signal

    """
    return numpy.real(cmdct(x, odd=odd)) * numpy.sqrt(2) 

Example 9

def ift(self):
      return MyImage(np.real(fft.ifft2(fft.fftshift(self.ft)))) 

Example 10

def ift(self):
        self.imgifft = MyImage(np.real(fft.ifft2(fft.fftshift(self.imgfft)))) 

Example 11

def get_polar_t(self):
        mag = self.get_magnitude()
        sizeimg = np.real(self.imgfft).shape
        
        pol = np.zeros(sizeimg)
        for x in range(sizeimg[0]):
            for y in range(sizeimg[1]):
                my = y - sizeimg[1] / 2
                mx = x - sizeimg[0] / 2
                if mx != 0:
                    phi = np.arctan(my / float(mx))
                else:
                    phi = 0
                r   = np.sqrt(mx**2 + my **2)
                
                ix = map_range(phi, -np.pi, np.pi, sizeimg[0], 0)
                iy = map_range(r, 0, sizeimg[0], 0, sizeimg[1])

                if ix >= 0 and ix < sizeimg[0] and iy >= 0 and iy < sizeimg[1]:
                    pol[x][y] =  mag.data[int(ix)][int(iy)]    
        pol = MyImage(pol)
        pol.limit(1)
        return pol 

Example 12

def operate(self, x):
        """
        Apply the separable filter to the signal vector *x*.
        """
        X = NP.fft.fftn(x, s=self.k_full)
        if NP.isrealobj(self.h) and NP.isrealobj(x):
            y = NP.real(NP.fft.ifftn(self.H * X))
        else:
            y = NP.fft.ifftn(self.H * X)

        if self.mode == 'full' or self.mode == 'circ':
            return y
        elif self.mode == 'valid':
            slice_list = []
            for i in range(self.ndim):
                if self.m[i]-1 == 0:
                    slice_list.append(slice(None, None, None))
                else:
                    slice_list.append(slice(self.m[i]-1, -(self.m[i]-1), None))
            return y[slice_list]
        else:
            assert(False) 

Example 13

def chain(mpas, astype=None):
    """Computes the tensor product of MPAs given in ``*args`` by adding more
    sites to the array.

    :param mpas: Iterable of MPAs in the order as they should appear in the
        chain
    :param astype: dtype of the returned MPA. If ``None``, use the type of the
        first MPA.
    :returns: MPA of length ``len(args[0]) + ... + len(args[-1])``

    .. todo:: Make this canonicalization aware
    .. todo:: Raise warning when casting complex to real dtype

    """
    mpas = iter(mpas)
    try:
        first = next(mpas)
    except StopIteration:
        raise ValueError('Argument `mpas` is an empty list')
    rest = (lt for mpa in mpas for lt in mpa.lt)
    if astype is None:
        astype = type(first)
    return astype(it.chain(first.lt, rest)) 

Example 14

def fftDf( df , part = "abs") :
   #Handle series or DataFrame
   if type(df) == pd.Series :
      df = pd.DataFrame(df)
      ise = True
   else :
      ise = False
   res = pd.DataFrame( index = np.fft.rfftfreq( df.index.size, d = dx( df ) ) )
   for col in df.columns :
      if part == "abs" :
         res["FFT_"+col] = np.abs( np.fft.rfft(df[col]) )  / (0.5*df.index.size)
      elif part == "real" :
         res["FFT_"+col] = np.real( np.fft.rfft(df[col]) ) / (0.5*df.index.size)
      elif part == "imag" :
         res["FFT_"+col] = np.imag( np.fft.rfft(df[col]) ) / (0.5*df.index.size)
   if ise :
      return res.iloc[:,0]
   else :
      return res 

Example 15

def derivFFT(df, n=1  ) :
   """ Deriv a signal trought FFT, warning, edge can be a bit noisy...
   indexList : channel to derive
   n : order of derivation
   """
   deriv = []
   for iSig in range(df.shape[1]) :
      fft = np.fft.fft( df.values[:,iSig] )   #FFT
      freq = np.fft.fftfreq( df.shape[0] , dx(df) )

      from copy import deepcopy
      fft0 = deepcopy(fft)
      if n>0 :
         fft *= (1j * 2*pi* freq[:])**n                    #Derivation in frequency domain
      else :
         fft[-n:] *= (1j * 2*pi* freq[-n:])**n
         fft[0:-n] = 0.

      tts = np.real(np.fft.ifft(fft))
      tts -= tts[0]
      deriv.append( tts )    #Inverse FFT

   return pd.DataFrame( data = np.transpose(deriv), index = df.index , columns = [ "DerivFFT("+ x +")" for x in df.columns ]  ) 

Example 16

def process(self, wave):
        wave.check_mono()
        if wave.sample_rate != self.sr:
            raise Exception("Wrong sample rate")                              
        n = int(np.ceil(2 * wave.num_frames / float(self.w_len)))
        m = (n + 1) * self.w_len / 2 
        swindow = self.make_signal_window(n)
        win_ratios = [self.window / swindow[t * self.w_len / 2 : 
            t * self.w_len / 2 + self.w_len] 
            for t in range(n)]
        wave = wave.zero_pad(0, int(m - wave.num_frames))
        wave = audio.Wave(signal.hilbert(wave), wave.sample_rate)        
        result = np.zeros((self.n_bins, n))
        
        for b in range(self.n_bins): 
            w = self.widths[b]
            wc = 1 / np.square(w + 1)
            filter = self.filters[b]
            band = fftfilt(filter, wave.zero_pad(0, int(2 * w))[:,0])
            band = band[int(w) : int(w + m), np.newaxis]    
            for t in range(n):
                frame = band[t * self.w_len / 2:
                             t * self.w_len / 2 + self.w_len,:] * win_ratios[t]
                result[b, t] =  wc * np.real(np.conj(np.dot(frame.conj().T, frame)))
        return audio.Spectrogram(result, self.sr, self.w_len, self.w_len / 2) 

Example 17

def test_psi(adjcube):
    """Tests retrieval of the wave functions and eigenvalues.
    """
    from pydft.bases.fourier import psi, O, H
    cell = adjcube
    V = QHO(cell)
    W = W4(cell)
    Ns = W.shape[1]
    Psi, epsilon = psi(V, W, cell, forceR=False)

    #Make sure that the eigenvalues are real.
    assert np.sum(np.imag(epsilon)) < 1e-13
    
    checkI = np.dot(Psi.conjugate().T, O(Psi, cell))
    assert abs(np.sum(np.diag(checkI))-Ns) < 1e-13 # Should be the identity
    assert np.abs(np.sum(checkI)-Ns) < 1e-13
    
    checkD = np.dot(Psi.conjugate().T, H(V, Psi, cell))
    diagsum = np.sum(np.diag(checkD))
    assert np.abs(np.sum(checkD)-diagsum) < 1e-12 # Should be diagonal

    # Should match the diagonal elements of previous matrix
    assert np.allclose(np.diag(checkD), epsilon) 

Example 18

def psi(V, W, cell, forceR=True):
    """Calculates the normalized wave functions using the basis coefficients.

    Args:
        V (pydft.potential.Potential): describing the potential for the
          particles.
        W (numpy.ndarray): wave function sample points.
        cell (pydft.geometry.Cell): describing the unit cell and sampling
          points.
        forceR (bool): forces the result to be real.
    """
    WN = Y(W, cell)
    mu = np.dot(WN.conjugate().T, H(V, WN, cell))
    epsilon, D = np.linalg.eig(mu)
    if forceR:
        epsilon = np.real(epsilon)
        
    return (np.dot(WN, D), epsilon) 

Example 19

def Idag(v=None, cell=None):
    """Computes the complex conjugate of the `I` operator for Fourier basis.

    Args:
        v (numpy.ndarray): if None, then return the matrix :math:`\mathbb{I^\dag}`,
          else, return :math:`\mathbb{I^\dag}\cdot v`.
        cell (pydft.geometry.Cell): that describes the unit cell and
          sampling points for real and reciprocal space.
    """
    #It turns out that for Fourier, the complex conjugate only differs by a -1
    #on the i (symmetric in R and G), so that we can return J instead.
    cell = get_cell(cell)
    def ifft(X):
        FB = np.fft.ifftn(np.reshape(X, cell.S, order='F'))
        return np.reshape(FB, X.shape, order='F')
    return matprod(ifft, v)*np.prod(cell.S) 

Example 20

def Jdag(v=None, cell=None):
    """Computes the complex conjugate of the `J` operator for Fourier basis.

    Args:
        v (numpy.ndarray): if None, then return the matrix :math:`\mathbb{J^\dag}`,
          else, return :math:`\mathbb{J^\dag}\cdot v`.
        cell (pydft.geometry.Cell): that describes the unit cell and
          sampling points for real and reciprocal space.
    """
    #It turns out that for Fourier, the complex conjugate only differs by a -1
    #on the i (symmetric in R and G), so that we can return I instead.
    cell = get_cell(cell)
    def fft(X):
        FF = np.fft.fftn(np.reshape(X, cell.S, order='F'))
        return np.reshape(FF, X.shape, order='F')
    return matprod(fft, v)/np.prod(cell.S) 

Example 21

def plotProfileXZplane(Ax,X,Z,Hx,Hz,Flag):

    FS = 20

    if Flag == 'Hp':
        Ax.streamplot(X,Z,Hx,Hz,color='b',linewidth=3.5,arrowsize=2)
        Ax.set_title('Primary Field',fontsize=FS+6)
    elif Flag == 'Hs_real':
        Ax.streamplot(X,Z,Hx,Hz,color='r',linewidth=3.5,arrowsize=2)
        Ax.set_title('Secondary Field (real)',fontsize=FS+6)
    elif Flag == 'Hs_imag':
        Ax.streamplot(X,Z,Hx,Hz,color='r',linewidth=3.5,arrowsize=2)
        Ax.set_title('Secondary Field (imaginary)',fontsize=FS+6)

    Ax.set_xbound(np.min(X),np.max(X))
    Ax.set_ybound(np.min(Z),np.max(Z))
    Ax.set_xlabel('X [m]',fontsize=FS+2)
    Ax.set_ylabel('Z [m]',fontsize=FS+2,labelpad=-10)
    Ax.tick_params(labelsize=FS-2) 

Example 22

def decode(ori_path, img_path, res_path, alpha):
    ori = cv2.imread(ori_path)
    img = cv2.imread(img_path)
    ori_f = np.fft.fft2(ori)
    img_f = np.fft.fft2(img)
    height, width = ori.shape[0], ori.shape[1]
    watermark = (ori_f - img_f) / alpha
    watermark = np.real(watermark)
    res = np.zeros(watermark.shape)
    random.seed(height + width)
    x = range(height / 2)
    y = range(width)
    random.shuffle(x)
    random.shuffle(y)
    for i in range(height / 2):
        for j in range(width):
            res[x[i]][y[j]] = watermark[i][j]
    cv2.imwrite(res_path, res, [int(cv2.IMWRITE_JPEG_QUALITY), 100]) 

Example 23

def genSpectra(time,dipole,signal):

    fw, frequency = pade(time,dipole)
    fw_sig, frequency = pade(time,signal,alternate=True)

    fw_re = np.real(fw)
    fw_im = np.imag(fw)
    fw_abs = fw_re**2 + fw_im**2

    #spectra = (fw_re*17.32)/(np.pi*field*damp_const)
    #spectra = (fw_re*17.32*514.220652)/(np.pi*field*damp_const)
    #numerator = np.imag((fw*np.conjugate(fw_sig)))
    numerator = np.imag(fw_abs*np.conjugate(fw_sig))
    #numerator = np.abs((fw*np.conjugate(fw_sig)))
    #numerator = np.abs(fw)
    denominator = np.real(np.conjugate(fw_sig)*fw_sig)
    #denominator = 1.0 
    spectra = ((4.0*27.21138602*2*frequency*np.pi*(numerator))/(3.0*137.036*denominator))
    spectra *= 1.0/100.0
    #plt.plot(frequency*27.2114,fourier)
    #plt.show()
    return frequency, spectra 

Example 24

def fft_convolve(X,Y, inv = 0):
    
    XF = np.fft.rfft2(X)
    YF = np.fft.rfft2(Y)
#    YF0 = np.copy(YF)
#    YF.imag = 0
#    XF.imag = 0
    if inv == 1:
 #       plt.imshow(np.real(YF)); plt.colorbar(); plt.show()
        YF = np.conj(YF)

    SF = XF*YF
    
    S = np.fft.irfft2(SF)
    n1,n2 = np.shape(S)

    S = np.roll(S,-n1/2+1,axis = 0)
    S = np.roll(S,-n2/2+1,axis = 1)

    return np.real(S) 

Example 25

def __init__(self,jet,kernels,k,x,y,pt,subpixel):
        self.jet = jet
        self.kernels = kernels
        self.k = k
        self.x = x
        self.y = y
        
        re = np.real(jet)
        im = np.imag(jet)
        
        self.mag = np.sqrt(re*re + im*im)
        self.phase = np.arctan2(re,im)
        
        if subpixel:
            d = np.array([[pt.X()-x],[pt.Y()-y]])
            comp = np.dot(self.k,d)
            self.phase -= comp.flatten()
            self.jet = self.mag*np.exp(1.0j*self.phase) 

Example 26

def calculate_katz_centrality(graph):
    """
    Compute the katz centrality for nodes.
    """
    # if not graph.is_directed():
    #    raise nx.NetworkXError( \
    #       "katz_centrality() not defined for undirected graphs.")
    print "\n\tCalculating Katz Centrality..."
    print "\tWarning: This might take a long time larger pedigrees."
    g = graph
    A = nx.adjacency_matrix(g)
    from scipy import linalg as LA
    max_eign = float(np.real(max(LA.eigvals(A.todense()))))
    print "\t-Max.Eigenvalue(A) ", round(max_eign, 3)
    kt = nx.katz_centrality(g, tol=1.0e-4, alpha=1/max_eign-0.01, beta=1.0, max_iter=999999)
    nx.set_node_attributes(g, 'katz', kt)
    katz_sorted = sorted(kt.items(), key=itemgetter(1), reverse=True)
    for key, value in katz_sorted[0:10]:
        print "\t   > ", key, round(value, 4)
    return g, kt 

Example 27

def pad_cpu2gpu(x, sz, offset=(0,0), dtype='real'):

    block_size = (16, 16 ,1)
    grid_size = (int(np.ceil(np.float32(sz[1])/block_size[1])),
                 int(np.ceil(np.float32(sz[0])/block_size[0])))

    sx = x.shape

    if x.__class__ == np.ndarray:
        x  = np.array(x).astype(np.float32)        
        x_gpu = cua.to_gpu(x)        
    elif x.__class__ == cua.GPUArray:       
        x_gpu = x

    if dtype == 'real':

        mod = cu.module_from_buffer(cubin)
        zeroPadKernel = mod.get_function("zeroPadKernel")

        x_padded_gpu = cua.zeros(tuple((int(sz[0]),int(sz[1]))), np.float32)
        
        zeroPadKernel(x_padded_gpu.gpudata, np.int32(sz[0]),     np.int32(sz[1]),
                             x_gpu.gpudata, np.int32(sx[0]),     np.int32(sx[1]),
                                            np.int32(offset[0]), np.int32(offset[1]),
                                            block=block_size, grid=grid_size)
    elif dtype == 'complex':

        mod = cu.module_from_buffer(cubin)
        #mod = SourceModule(open('gputools.cu').read(), keep=True)
        zeroPadComplexKernel = mod.get_function("zeroPadComplexKernel")

        x_padded_gpu = cua.zeros(tuple((int(sz[0]),int(sz[1]))), np.complex64)
        
        zeroPadComplexKernel(x_padded_gpu.gpudata, np.int32(sz[0]),     np.int32(sz[1]),
                                    x_gpu.gpudata, np.int32(sx[0]),     np.int32(sx[1]),
                                                   np.int32(offset[0]), np.int32(offset[1]),
                                                   block=block_size, grid=grid_size)

    return x_padded_gpu 

Example 28

def cnvinv_gradfun(self, z, sz, y_gpu, alpha=0., beta=0.):
        """
        Computes gradient used for 'lbfgsb' mode of deconv method.
        See deconv for details.
        """
        
        if z.__class__ == np.ndarray:
            z = np.array(np.reshape(z,sz)).astype(np.float32)
            z_gpu = cua.to_gpu(z)            

        grad_gpu =  self.cnvtp(self.res_gpu)
        
        # Thikonov regularization
        # alpha > 0: Thikonov on the gradient of z
        if alpha > 0:
            grad_gpu += alpha * self.lz_gpu

        # beta  > 0: Thikonov on z
        if beta > 0:
            grad_gpu += beta * z_gpu

        grad = -np.real(grad_gpu.get())
        grad = grad.flatten()
        return grad.astype(np.float64) 

Example 29

def psfScale(D, wavelength, pixSize):
  """
  Return the PSF scale appropriate for the required pixel size, wavelength and telescope diameter
  The aperture is padded by this amount; resultant pix scale is lambda/D/psf_scale, so for instance full frame 256 pix
  for 3.5 m at 532 nm is 256*5.32e-7/3.5/3 = 2.67 arcsec for psf_scale = 3
  
  Args:
      D (real): telescope diameter in m
      wavelength (real): wavelength in Angstrom
      pixSize (real): pixel size in arcsec
  
  Returns:
      real: psf scale
  """
  DInCm = D * 100.0
  wavelengthInCm = wavelength * 1e-8
  return 206265.0 * wavelengthInCm / (DInCm * pixSize) 

Example 30

def fwd_filter(img, S):
    img_w, img_h, ch = img.shape
    F = pad2(img)
    F.data[F.mask] = 0.  # make sure its zero-filled!

    # Forward transform
    specF = np.fft.fft2(F.data.astype(float), axes=(0, 1))
    specN = np.fft.fft2(1. - F.mask.astype(float), axes=(0, 1))
    specS = np.fft.fft2(S[::-1, ::-1])
    out = np.real(np.fft.ifft2(specF * specS[:, :, np.newaxis], axes=(0, 1)))
    norm = np.real(np.fft.ifft2(specN * specS[:, :, np.newaxis], axes=(0, 1)))
    eps = 1e-15
    norm = np.maximum(norm, eps)
    out /= norm
    out = out[-img_w:, -img_h:]
    out[img.mask] = 0.
    return np.ma.MaskedArray(data=out, mask=img.mask) 

Example 31

def kernel_impute(img, S):
    F = pad2(img)
    F.data[F.mask] = 0.  # make sure its zero-filled!
    img_w, img_h, img_ch = img.shape
    Q = S
    specF = np.fft.fft2(F.data.astype(float), axes=(0, 1))
    specN = np.fft.fft2(1. - F.mask.astype(float), axes=(0, 1))
    specQ = np.fft.fft2(Q[::-1, ::-1])
    numer = np.real(np.fft.ifft2(specF * specQ[:, :, np.newaxis], axes=(0, 1)))
    denom = np.real(np.fft.ifft2(specN * specQ[:, :, np.newaxis], axes=(0, 1)))
    eps = 1e-15
    fill = numer/(denom+eps)
    fill = fill[-img_w:, -img_h:]

    image = img.data.copy()

    # img = img.copy()
    image[img.mask] = fill[img.mask]
    mask = np.zeros_like(img.mask, dtype=bool)
    return np.ma.MaskedArray(data=image, mask=mask) 

Example 32

def get_neg_log_post(Phi, sigma_J_list, ROI_list, G, MMT, q, Sigma_E,  GL,
                     nu, V, prior_on = False):
    eps = 1E-13
    p = Phi.shape[0]
    n_ROI = len(sigma_J_list)
    Qu = Phi.dot(Phi.T)
    G_Sigma_G = np.zeros(MMT.shape)
    for i in range(n_ROI):
        G_Sigma_G += sigma_J_list[i]**2 * np.dot(G[:,ROI_list[i]], G[:,ROI_list[i]].T)
    cov = Sigma_E + G_Sigma_G + GL.dot(Qu).dot(GL.T) 
    inv_cov = np.linalg.inv(cov)    
    eigs = np.real(np.linalg.eigvals(cov)) + eps
    log_det_cov = np.sum(np.log(eigs))  
    result = q*log_det_cov + np.trace(MMT.dot(inv_cov))
    if prior_on:
        inv_Q = np.linalg.inv(Qu)
        #det_Q = np.linalg.det(Qu)
        log_det_Q = np.sum(np.log(np.diag(Phi)**2))
        result =  result + np.float(nu+p+1)*log_det_Q+ np.trace(V.dot(inv_Q))
    return result

#==============================================================================
# update both Qu and Sigma_J, gradient of Qu and Sigma J 

Example 33

def get_neg_log_post(Phi, sigma_J_list, ROI_list, G, MMT, q, Sigma_E,  GL,
                     nu, V, prior_on = False):
    eps = 1E-13
    p = Phi.shape[0]
    n_ROI = len(sigma_J_list)
    Qu = Phi.dot(Phi.T)
    G_Sigma_G = np.zeros(MMT.shape)
    for i in range(n_ROI):
        G_Sigma_G += sigma_J_list[i]**2 * np.dot(G[:,ROI_list[i]], G[:,ROI_list[i]].T)
    cov = Sigma_E + G_Sigma_G + GL.dot(Qu).dot(GL.T) 
    inv_cov = np.linalg.inv(cov)    
    eigs = np.real(np.linalg.eigvals(cov)) + eps
    log_det_cov = np.sum(np.log(eigs))  
    result = q*log_det_cov + np.trace(MMT.dot(inv_cov))
    if prior_on:
        inv_Q = np.linalg.inv(Qu)
        #det_Q = np.linalg.det(Qu)
        log_det_Q = np.sum(np.log(np.diag(Phi)**2))
        result =  result + np.float(nu+p+1)*log_det_Q+ np.trace(V.dot(inv_Q))
    return result

#==============================================================================
# update both Qu and Sigma_J, gradient of Qu and Sigma J 

Example 34

def __init__(self, qubit_names, quad="real"):
        super(PulseCalibration, self).__init__()
        self.qubit_names = qubit_names if isinstance(qubit_names, list) else [qubit_names]
        self.qubit     = [QubitFactory(qubit_name) for qubit_name in qubit_names] if isinstance(qubit_names, list) else QubitFactory(qubit_names)
        self.filename   = 'None'
        self.exp        = None
        self.axis_descriptor = None
        self.cw_mode    = False
        self.saved_settings = config.load_meas_file(config.meas_file)
        self.settings = deepcopy(self.saved_settings) #make a copy for used during calibration
        self.quad = quad
        if quad == "real":
            self.quad_fun = np.real
        elif quad == "imag":
            self.quad_fun = np.imag
        elif quad == "amp":
            self.quad_fun = np.abs
        elif quad == "phase":
            self.quad_fun = np.angle
        else:
            raise ValueError('Quadrature to calibrate must be one of ("real", "imag", "amp", "phase").')
        self.plot       = self.init_plot() 

Example 35

def find_null_offset(xpts, powers, default=0.0):
    """Finds the offset corresponding to the minimum power using a fit to the measured data"""
    def model(x, a, b, c):
        return a*(x - b)**2 + c
    powers = np.power(10, powers/10.)
    min_idx = np.argmin(powers)
    try:
        fit = curve_fit(model, xpts, powers, p0=[1, xpts[min_idx], powers[min_idx]])
    except RuntimeError:
        logger.warning("Mixer null offset fit failed.")
        return default, np.zeros(len(powers))
    best_offset = np.real(fit[0][1])
    best_offset = np.minimum(best_offset, xpts[-1])
    best_offset = np.maximum(best_offset, xpts[0])
    xpts_fine = np.linspace(xpts[0],xpts[-1],101)
    fit_pts = np.array([np.real(model(x, *fit[0])) for x in xpts_fine])
    if min(fit_pts)<0: fit_pts-=min(fit_pts)-1e-10 #prevent log of a negative number
    return best_offset, xpts_fine, 10*np.log10(fit_pts) 

Example 36

def make_layout(self):
        self.lay = QtWidgets.QHBoxLayout()
        self.lay.setContentsMargins(0, 0, 0, 0)
        self.real = FloatSpinBox(label=self.labeltext,
                                 min=self.minimum,
                                 max=self.maximum,
                                 increment=self.singleStep,
                                 log_increment=self.log_increment,
                                 halflife_seconds=self.halflife_seconds,
                                 decimals=self.decimals)
        self.imag = FloatSpinBox(label=self.labeltext,
                                 min=self.minimum,
                                 max=self.maximum,
                                 increment=self.singleStep,
                                 log_increment=self.log_increment,
                                 halflife_seconds=self.halflife_seconds,
                                 decimals=self.decimals)
        self.real.value_changed.connect(self.value_changed)
        self.lay.addWidget(self.real)
        self.label = QtWidgets.QLabel(" + j")
        self.lay.addWidget(self.label)
        self.imag.value_changed.connect(self.value_changed)
        self.lay.addWidget(self.imag)
        self.setLayout(self.lay)
        self.setFocusPolicy(QtCore.Qt.ClickFocus) 

Example 37

def save_image(imager, grid_data, grid_norm, output_file):
    """Makes an image from gridded visibilities and saves it to a FITS file.

    Args:
        imager (oskar.Imager):          Handle to configured imager.
        grid_data (numpy.ndarray):      Final visibility grid.
        grid_norm (float):              Grid normalisation to apply.
        output_file (str):              Name of output FITS file to write.
    """
    # Make the image (take the FFT, normalise, and apply grid correction).
    imager.finalise_plane(grid_data, grid_norm)
    grid_data = numpy.real(grid_data)

    # Trim the image if required.
    border = (imager.plane_size - imager.image_size) // 2
    if border > 0:
        end = border + imager.image_size
        grid_data = grid_data[border:end, border:end]

    # Write the FITS file.
    hdr = fits.header.Header()
    fits.writeto(output_file, grid_data, hdr, clobber=True) 

Example 38

def show_image(im: Image, fig=None, title: str = '', pol=0, chan=0, cm='rainbow'):
    """ Show an Image with coordinates using matplotlib

    :param im:
    :param fig:
    :param title:
    :return:
    """
    import matplotlib.pyplot as plt

    assert isinstance(im, Image)
    if not fig:
        fig = plt.figure()
    plt.clf()
    fig.add_subplot(111, projection=im.wcs.sub(['longitude', 'latitude']))
    if len(im.data.shape) == 4:
        plt.imshow(numpy.real(im.data[chan, pol, :, :]), origin='lower', cmap=cm)
    elif len(im.data.shape) == 2:
        plt.imshow(numpy.real(im.data[:, :]), origin='lower', cmap=cm)
    plt.xlabel('RA---SIN')
    plt.ylabel('DEC--SIN')
    plt.title(title)
    plt.colorbar()
    return fig 

Example 39

def convolve_convolve_scalestack(scalestack, img):
    """Convolve img by the specified scalestack, returning the resulting stack

    :param scalestack: stack containing the scales
    :param img: Image to be convolved
    :return: Twice convolved image [nscales, nscales, nx, ny]
    """

    nscales, nx, ny = scalestack.shape
    convolved_shape = [nscales, nscales, nx, ny]
    convolved = numpy.zeros(convolved_shape)
    ximg = numpy.fft.fftshift(numpy.fft.fft2(numpy.fft.fftshift(img)))

    xscaleshape = [nscales, nx, ny]
    xscale = numpy.zeros(xscaleshape, dtype='complex')
    for s in range(nscales):
        xscale[s] = numpy.fft.fftshift(numpy.fft.fft2(numpy.fft.fftshift(scalestack[s, ...])))

    for s in range(nscales):
        for p in range(nscales):
            xmult = ximg * xscale[p] * numpy.conjugate(xscale[s])
            convolved[s, p, ...] = numpy.real(numpy.fft.ifftshift(numpy.fft.ifft2(numpy.fft.ifftshift(xmult))))
    return convolved 

Example 40

def plot_waveforms(waveforms, figTitle=''):
    channels = waveforms.keys()
    # plot
    plots = []
    for (ct, chan) in enumerate(channels):
        fig = bk.figure(title=figTitle + repr(chan),
                        plot_width=800,
                        plot_height=350,
                        y_range=[-1.05, 1.05],
                        x_axis_label=u'Time (?s)')
        fig.background_fill_color = config.plotBackground
        if config.gridColor:
            fig.xgrid.grid_line_color = config.gridColor
            fig.ygrid.grid_line_color = config.gridColor
        waveformToPlot = waveforms[chan]
        xpts = np.linspace(0, len(waveformToPlot) / chan.phys_chan.sampling_rate
                           / 1e-6, len(waveformToPlot))
        fig.line(xpts, np.real(waveformToPlot), color='red')
        fig.line(xpts, np.imag(waveformToPlot), color='blue')
        plots.append(fig)
    bk.show(column(*plots)) 

Example 41

def merge_waveform(n, chAB, chAm1, chAm2, chBm1, chBm2):
    '''
    Builds packed I and Q waveforms from the nth mini LL, merging in marker data.
    '''
    wfAB = np.array([], dtype=np.complex)
    for entry in chAB['linkList'][n % len(chAB['linkList'])]:
        if not entry.isTimeAmp:
            wfAB = np.append(wfAB, chAB['wfLib'][entry.key])
        else:
            wfAB = np.append(wfAB, chAB['wfLib'][entry.key][0] *
                             np.ones(entry.length * entry.repeat))

    wfAm1 = marker_waveform(chAm1['linkList'][n % len(chAm1['linkList'])],
                            chAm1['wfLib'])
    wfAm2 = marker_waveform(chAm2['linkList'][n % len(chAm2['linkList'])],
                            chAm2['wfLib'])
    wfBm1 = marker_waveform(chBm1['linkList'][n % len(chBm1['linkList'])],
                            chBm1['wfLib'])
    wfBm2 = marker_waveform(chBm2['linkList'][n % len(chBm2['linkList'])],
                            chBm2['wfLib'])

    wfA = pack_waveform(np.real(wfAB), wfAm1, wfAm2)
    wfB = pack_waveform(np.imag(wfAB), wfBm1, wfBm2)

    return wfA, wfB 

Example 42

def tune_everything(x0squared, C, T, gmin, gmax):
  # First tune based on dynamic range    
  if C==0:
    dr=gmax/gmin
    mustar=((np.sqrt(dr)-1)/(np.sqrt(dr)+1))**2
    alpha_star = (1+np.sqrt(mustar))**2/gmax
    
    return alpha_star,mustar

  dist_to_opt = x0squared
  grad_var = C
  max_curv = gmax
  min_curv = gmin
  const_fact = dist_to_opt * min_curv**2 / 2 / grad_var
  coef = [-1, 3, -(3 + const_fact), 1]
  roots = np.roots(coef)
  roots = roots[np.real(roots) > 0]
  roots = roots[np.real(roots) < 1]
  root = roots[np.argmin(np.imag(roots) ) ]

  assert root > 0 and root < 1 and np.absolute(root.imag) < 1e-6

  dr = max_curv / min_curv
  assert max_curv >= min_curv
  mu = max( ( (np.sqrt(dr) - 1) / (np.sqrt(dr) + 1) )**2, root**2)

  lr_min = (1 - np.sqrt(mu) )**2 / min_curv
  lr_max = (1 + np.sqrt(mu) )**2 / max_curv

  alpha_star = lr_min
  mustar = mu

  return alpha_star, mustar 

Example 43

def reflection_from_matrix(matrix):
    """Return mirror plane point and normal vector from reflection matrix.

    >>> v0 = numpy.random.random(3) - 0.5
    >>> v1 = numpy.random.random(3) - 0.5
    >>> M0 = reflection_matrix(v0, v1)
    >>> point, normal = reflection_from_matrix(M0)
    >>> M1 = reflection_matrix(point, normal)
    >>> is_same_transform(M0, M1)
    True

    """
    M = numpy.array(matrix, dtype=numpy.float64, copy=False)
    # normal: unit eigenvector corresponding to eigenvalue -1
    l, V = numpy.linalg.eig(M[:3, :3])
    i = numpy.where(abs(numpy.real(l) + 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue -1")
    normal = numpy.real(V[:, i[0]]).squeeze()
    # point: any unit eigenvector corresponding to eigenvalue 1
    l, V = numpy.linalg.eig(M)
    i = numpy.where(abs(numpy.real(l) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    point = numpy.real(V[:, i[-1]]).squeeze()
    point /= point[3]
    return point, normal 

Example 44

def rotation_from_matrix(matrix):
    """Return rotation angle and axis from rotation matrix.

    >>> angle = (random.random() - 0.5) * (2*math.pi)
    >>> direc = numpy.random.random(3) - 0.5
    >>> point = numpy.random.random(3) - 0.5
    >>> R0 = rotation_matrix(angle, direc, point)
    >>> angle, direc, point = rotation_from_matrix(R0)
    >>> R1 = rotation_matrix(angle, direc, point)
    >>> is_same_transform(R0, R1)
    True

    """
    R = numpy.array(matrix, dtype=numpy.float64, copy=False)
    R33 = R[:3, :3]
    # direction: unit eigenvector of R33 corresponding to eigenvalue of 1
    l, W = numpy.linalg.eig(R33.T)
    i = numpy.where(abs(numpy.real(l) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    direction = numpy.real(W[:, i[-1]]).squeeze()
    # point: unit eigenvector of R33 corresponding to eigenvalue of 1
    l, Q = numpy.linalg.eig(R)
    i = numpy.where(abs(numpy.real(l) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    point = numpy.real(Q[:, i[-1]]).squeeze()
    point /= point[3]
    # rotation angle depending on direction
    cosa = (numpy.trace(R33) - 1.0) / 2.0
    if abs(direction[2]) > 1e-8:
        sina = (R[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2]
    elif abs(direction[1]) > 1e-8:
        sina = (R[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1]
    else:
        sina = (R[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0]
    angle = math.atan2(sina, cosa)
    return angle, direction, point 

Example 45

def scale_from_matrix(matrix):
    """Return scaling factor, origin and direction from scaling matrix.

    >>> factor = random.random() * 10 - 5
    >>> origin = numpy.random.random(3) - 0.5
    >>> direct = numpy.random.random(3) - 0.5
    >>> S0 = scale_matrix(factor, origin)
    >>> factor, origin, direction = scale_from_matrix(S0)
    >>> S1 = scale_matrix(factor, origin, direction)
    >>> is_same_transform(S0, S1)
    True
    >>> S0 = scale_matrix(factor, origin, direct)
    >>> factor, origin, direction = scale_from_matrix(S0)
    >>> S1 = scale_matrix(factor, origin, direction)
    >>> is_same_transform(S0, S1)
    True

    """
    M = numpy.array(matrix, dtype=numpy.float64, copy=False)
    M33 = M[:3, :3]
    factor = numpy.trace(M33) - 2.0
    try:
        # direction: unit eigenvector corresponding to eigenvalue factor
        l, V = numpy.linalg.eig(M33)
        i = numpy.where(abs(numpy.real(l) - factor) < 1e-8)[0][0]
        direction = numpy.real(V[:, i]).squeeze()
        direction /= vector_norm(direction)
    except IndexError:
        # uniform scaling
        factor = (factor + 2.0) / 3.0
        direction = None
    # origin: any eigenvector corresponding to eigenvalue 1
    l, V = numpy.linalg.eig(M)
    i = numpy.where(abs(numpy.real(l) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no eigenvector corresponding to eigenvalue 1")
    origin = numpy.real(V[:, i[-1]]).squeeze()
    origin /= origin[3]
    return factor, origin, direction 

Example 46

def reflection_from_matrix(matrix):
    """Return mirror plane point and normal vector from reflection matrix.

    >>> v0 = numpy.random.random(3) - 0.5
    >>> v1 = numpy.random.random(3) - 0.5
    >>> M0 = reflection_matrix(v0, v1)
    >>> point, normal = reflection_from_matrix(M0)
    >>> M1 = reflection_matrix(point, normal)
    >>> is_same_transform(M0, M1)
    True

    """
    M = numpy.array(matrix, dtype=numpy.float64, copy=False)
    # normal: unit eigenvector corresponding to eigenvalue -1
    w, V = numpy.linalg.eig(M[:3, :3])
    i = numpy.where(abs(numpy.real(w) + 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue -1")
    normal = numpy.real(V[:, i[0]]).squeeze()
    # point: any unit eigenvector corresponding to eigenvalue 1
    w, V = numpy.linalg.eig(M)
    i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    point = numpy.real(V[:, i[-1]]).squeeze()
    point /= point[3]
    return point, normal 

Example 47

def rotation_from_matrix(matrix):
    """Return rotation angle and axis from rotation matrix.

    >>> angle = (random.random() - 0.5) * (2*math.pi)
    >>> direc = numpy.random.random(3) - 0.5
    >>> point = numpy.random.random(3) - 0.5
    >>> R0 = rotation_matrix(angle, direc, point)
    >>> angle, direc, point = rotation_from_matrix(R0)
    >>> R1 = rotation_matrix(angle, direc, point)
    >>> is_same_transform(R0, R1)
    True

    """
    R = numpy.array(matrix, dtype=numpy.float64, copy=False)
    R33 = R[:3, :3]
    # direction: unit eigenvector of R33 corresponding to eigenvalue of 1
    w, W = numpy.linalg.eig(R33.T)
    i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    direction = numpy.real(W[:, i[-1]]).squeeze()
    # point: unit eigenvector of R33 corresponding to eigenvalue of 1
    w, Q = numpy.linalg.eig(R)
    i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    point = numpy.real(Q[:, i[-1]]).squeeze()
    point /= point[3]
    # rotation angle depending on direction
    cosa = (numpy.trace(R33) - 1.0) / 2.0
    if abs(direction[2]) > 1e-8:
        sina = (R[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2]
    elif abs(direction[1]) > 1e-8:
        sina = (R[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1]
    else:
        sina = (R[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0]
    angle = math.atan2(sina, cosa)
    return angle, direction, point 

Example 48

def scale_from_matrix(matrix):
    """Return scaling factor, origin and direction from scaling matrix.

    >>> factor = random.random() * 10 - 5
    >>> origin = numpy.random.random(3) - 0.5
    >>> direct = numpy.random.random(3) - 0.5
    >>> S0 = scale_matrix(factor, origin)
    >>> factor, origin, direction = scale_from_matrix(S0)
    >>> S1 = scale_matrix(factor, origin, direction)
    >>> is_same_transform(S0, S1)
    True
    >>> S0 = scale_matrix(factor, origin, direct)
    >>> factor, origin, direction = scale_from_matrix(S0)
    >>> S1 = scale_matrix(factor, origin, direction)
    >>> is_same_transform(S0, S1)
    True

    """
    M = numpy.array(matrix, dtype=numpy.float64, copy=False)
    M33 = M[:3, :3]
    factor = numpy.trace(M33) - 2.0
    try:
        # direction: unit eigenvector corresponding to eigenvalue factor
        w, V = numpy.linalg.eig(M33)
        i = numpy.where(abs(numpy.real(w) - factor) < 1e-8)[0][0]
        direction = numpy.real(V[:, i]).squeeze()
        direction /= vector_norm(direction)
    except IndexError:
        # uniform scaling
        factor = (factor + 2.0) / 3.0
        direction = None
    # origin: any eigenvector corresponding to eigenvalue 1
    w, V = numpy.linalg.eig(M)
    i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no eigenvector corresponding to eigenvalue 1")
    origin = numpy.real(V[:, i[-1]]).squeeze()
    origin /= origin[3]
    return factor, origin, direction 

Example 49

def rotation_from_matrix(matrix):
    """Return rotation angle and axis from rotation matrix.

    >>> angle = (random.random() - 0.5) * (2*math.pi)
    >>> direc = numpy.random.random(3) - 0.5
    >>> point = numpy.random.random(3) - 0.5
    >>> R0 = rotation_matrix(angle, direc, point)
    >>> angle, direc, point = rotation_from_matrix(R0)
    >>> R1 = rotation_matrix(angle, direc, point)
    >>> is_same_transform(R0, R1)
    True

    """
    R = numpy.array(matrix, dtype=numpy.float64, copy=False)
    R33 = R[:3, :3]
    # direction: unit eigenvector of R33 corresponding to eigenvalue of 1
    w, W = numpy.linalg.eig(R33.T)
    i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    direction = numpy.real(W[:, i[-1]]).squeeze()
    # point: unit eigenvector of R33 corresponding to eigenvalue of 1
    w, Q = numpy.linalg.eig(R)
    i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    point = numpy.real(Q[:, i[-1]]).squeeze()
    point /= point[3]
    # rotation angle depending on direction
    cosa = (numpy.trace(R33) - 1.0) / 2.0
    if abs(direction[2]) > 1e-8:
        sina = (R[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2]
    elif abs(direction[1]) > 1e-8:
        sina = (R[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1]
    else:
        sina = (R[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0]
    angle = math.atan2(sina, cosa)
    return angle, direction, point

# Function to translate handshape coding to degrees of rotation, adduction, flexion 

Example 50

def rotation_from_matrix(matrix):
    """Return rotation angle and axis from rotation matrix.

    >>> angle = (random.random() - 0.5) * (2*math.pi)
    >>> direc = numpy.random.random(3) - 0.5
    >>> point = numpy.random.random(3) - 0.5
    >>> R0 = rotation_matrix(angle, direc, point)
    >>> angle, direc, point = rotation_from_matrix(R0)
    >>> R1 = rotation_matrix(angle, direc, point)
    >>> is_same_transform(R0, R1)
    True

    """
    R = numpy.array(matrix, dtype=numpy.float64, copy=False)
    R33 = R[:3, :3]
    # direction: unit eigenvector of R33 corresponding to eigenvalue of 1
    w, W = numpy.linalg.eig(R33.T)
    i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    direction = numpy.real(W[:, i[-1]]).squeeze()
    # point: unit eigenvector of R33 corresponding to eigenvalue of 1
    w, Q = numpy.linalg.eig(R)
    i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
    if not len(i):
        raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
    point = numpy.real(Q[:, i[-1]]).squeeze()
    point /= point[3]
    # rotation angle depending on direction
    cosa = (numpy.trace(R33) - 1.0) / 2.0
    if abs(direction[2]) > 1e-8:
        sina = (R[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2]
    elif abs(direction[1]) > 1e-8:
        sina = (R[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1]
    else:
        sina = (R[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0]
    angle = math.atan2(sina, cosa)
    return angle, direction, point

# Function to translate handshape coding to degrees of rotation, adduction, flexion 
点赞

发表评论

电子邮件地址不会被公开。 必填项已用*标注