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])
python类real()的实例源码
transformations.py 文件源码
项目:Neural-Networks-for-Inverse-Kinematics
作者: paramrajpura
项目源码
文件源码
阅读 47
收藏 0
点赞 0
评论 0
def quaternion_real(quaternion):
"""Return real part of quaternion.
>>> quaternion_real([3, 0, 1, 2])
3.0
"""
return float(quaternion[0])
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
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
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
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)
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
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
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)
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))
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
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 ] )
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)
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)
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)
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)
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)
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)
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])
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
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)
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)
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
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)
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)
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)
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)
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
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
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()