def flow2parallax(u,v,q):
"""
Given the flow fields (after correction!) and the epipole,
return:
- The normalized parallax (HxW array)
- The vectors pointing to the epipoles (HxWx2 array)
- The distances of all points to the epipole (HxW array)
"""
h,w = u.shape
y,x = np.mgrid[:h,:w]
u_f = q[0] - x
v_f = q[1] - y
dists = np.sqrt(u_f**2 + v_f**2)
u_f_n = u_f / np.maximum(dists,1e-3)
v_f_n = v_f / np.maximum(dists,1e-3)
parallax = u * u_f_n + v * v_f_n
return parallax, np.dstack((u_f_n, v_f_n)), dists
python类mgrid()的实例源码
def create_test_dataset(image_shape, n, circle_radius, donut_radius):
img = np.zeros((image_shape[0], image_shape[1]))
y_pixels = np.arange(0, image_shape[0], 1)
x_pixels = np.arange(0, image_shape[1], 1)
cell_y_coords = np.random.choice(y_pixels, n, replace=False)
cell_x_coords = np.random.choice(x_pixels, n, replace=False)
for x, y in zip(cell_x_coords, cell_y_coords):
xx, yy = np.mgrid[:512, :512] # create mesh grid of image dimensions
circle = (xx - x) ** 2 + (yy - y) ** 2 # apply circle formula
donut = np.logical_and(circle < (circle_radius+donut_radius),
circle > (circle_radius-5)) # donuts are thresholded circles
thresholded_circle = circle < circle_radius
img[np.where(thresholded_circle)] = 1
img[np.where(donut)] = 2
return img
def test_square_grid():
X = np.mgrid[0:16, 0:16]
X = X.reshape((len(X), -1)).T
name = 'square'
D, Q = test_toy_embedding(X, 32, 2, name, palette='hls')
def plot_mat_on_data(mat, sample):
plt.figure()
plot_data_embedded(X, palette='w')
alpha = np.maximum(mat[sample], 0) / mat[sample].max()
plot_data_embedded(X, palette='#FF0000', alpha=alpha)
pdf_file_name = '{}{}_plot_{}_on_data_{}{}'
plot_mat_on_data(D, 7 * 16 + 7)
plt.savefig(pdf_file_name.format(dir_name, name, 'D', 'middle', '.pdf'))
plot_mat_on_data(Q, 7 * 16 + 7)
plt.savefig(pdf_file_name.format(dir_name, name, 'Q', 'middle', '.pdf'))
# for s in range(len(X)):
# plot_mat_on_data(Q, s)
# plt.savefig(pdf_file_name.format(dir_name, name, 'Q', s, '.png'))
# plt.close()
def plot(self, ax, idx1, idx2, range1, range2, n=100):
assert len(range1) == len(range2) == 2 and idx1 != idx2
x, y = np.mgrid[range1[0]:range1[1]:(n+0j), range2[0]:range2[1]:(n+0j)]
if isinstance(self.action_space, ContinuousSpace):
points_B_Doa = np.zeros((n*n, self.obsfeat_space.storage_size + self.action_space.storage_size))
points_B_Doa[:,idx1] = x.ravel()
points_B_Doa[:,idx2] = y.ravel()
obsfeat_B_Df, a_B_Da = points_B_Doa[:,:self.obsfeat_space.storage_size], points_B_Doa[:,self.obsfeat_space.storage_size:]
assert a_B_Da.shape[1] == self.action_space.storage_size
t_B = np.zeros(a_B_Da.shape[0]) # XXX make customizable
z = self.compute_reward(obsfeat_B_Df, a_B_Da, t_B).reshape(x.shape)
else:
obsfeat_B_Df = np.zeros((n*n, self.obsfeat_space.storage_size))
obsfeat_B_Df[:,idx1] = x.ravel()
obsfeat_B_Df[:,idx2] = y.ravel()
a_B_Da = np.zeros((obsfeat_B_Df.shape[0], 1), dtype=np.int32) # XXX make customizable
t_B = np.zeros(a_B_Da.shape[0]) # XXX make customizable
z = self.compute_reward(obsfeat_B_Df, a_B_Da, t_B).reshape(x.shape)
ax.pcolormesh(x, y, z, cmap='viridis')
ax.contour(x, y, z, levels=np.log(np.linspace(2., 3., 10)))
# ax.contourf(x, y, z, levels=[np.log(2.), np.log(2.)+.5], alpha=.5) # high-reward region is highlighted
def calculate_scalar_matrix(values_a, values_b):
"""
convenience function wrapper of py:function:`calculate_scalar_product_matrix` for the case of scalar elements.
:param values_a:
:param values_b:
:return:
"""
return calculate_scalar_product_matrix(np.multiply,
sanitize_input(values_a, Number),
sanitize_input(values_b, Number))
# i, j = np.mgrid[0:values_a.shape[0], 0:values_b.shape[0]]
# vals_i = values_a[i]
# vals_j = values_b[j]
# return np.multiply(vals_i, vals_j)
def gabor_2d(M, N, sigma, theta, xi, slant=1.0, offset=0, fft_shift=None):
gab = np.zeros((M, N), np.complex64)
R = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]], np.float32)
R_inv = np.array([[np.cos(theta), np.sin(theta)], [-np.sin(theta), np.cos(theta)]], np.float32)
D = np.array([[1, 0], [0, slant * slant]])
curv = np.dot(R, np.dot(D, R_inv)) / ( 2 * sigma * sigma)
for ex in [-2, -1, 0, 1, 2]:
for ey in [-2, -1, 0, 1, 2]:
[xx, yy] = np.mgrid[offset + ex * M:offset + M + ex * M, offset + ey * N:offset + N + ey * N]
arg = -(curv[0, 0] * np.multiply(xx, xx) + (curv[0, 1] + curv[1, 0]) * np.multiply(xx, yy) + curv[
1, 1] * np.multiply(yy, yy)) + 1.j * (xx * xi * np.cos(theta) + yy * xi * np.sin(theta))
gab = gab + np.exp(arg)
norm_factor = (2 * 3.1415 * sigma * sigma / slant)
gab = gab / norm_factor
if (fft_shift):
gab = np.fft.fftshift(gab, axes=(0, 1))
return gab
def test():
import PIL.Image
y, x = np.mgrid[0:256, 0:256]
z = np.ones((256,256)) * 128
img0 = np.dstack((x, y, z)).astype(np.uint8)
img1 = y.astype(np.uint8)
img2 = np.arange(256, dtype=np.uint8)
img3 = PIL.Image.open("pics/RGB.png")
img3 = np.array(img3)[:,:,0:3]
img4 = PIL.Image.open("pics/banff.jpg")
img4 = np.array(img4)[:,:,0:3]
img5, _ = (np.mgrid[0:1242, 0:1276] / 1242. * 255.).astype(np.uint8)
img6, _ = (np.mgrid[0:1007, 0:12] / 1007. * 255.).astype(np.uint8)
for i in (1, 2, 4, 8):
write_tiff("Test0_" + str(i) + ".TIF", img0, bit_depth=i)
write_tiff("Test1_" + str(i) + ".TIF", img1, bit_depth=i)
write_tiff("Test2_" + str(i) + ".TIF", img2, bit_depth=i)
write_tiff("Test3_" + str(i) + ".TIF", img3, bit_depth=i)
write_tiff("Test4_" + str(i) + ".TIF", img4, bit_depth=i)
write_tiff("Test5_" + str(i) + ".TIF", img5, bit_depth=i)
write_tiff("Test6_" + str(i) + ".TIF", img6, bit_depth=i)
def gauss_kernel(size, sigma=None, size_y=None, sigma_y=None):
"""
Generates a 2D Gaussian kernel as a numpy array
Args:
size (int): 1/2 the width of the kernel; total width := 2*size+1
sigma (float): spread of the gaussian in the width direction
size_y (int): 1/2 the height of the kernel; defaults to size
sigma_y (float): spread of the gaussian in the height direction; defaults to sigma
Returns:
numpy array: normalized 2D gaussian array
"""
size = int(size)
if not size_y:
size_y = size
else:
size_y = int(size_y)
if not sigma:
sigma = 0.5 * size + .1
if not sigma_y:
sigma_y = sigma
x, y = np.mgrid[-size:size+1, -size_y:size_y+1]
g = np.exp(-0.5 * (x ** 2 / sigma ** 2 + y ** 2 / sigma_y ** 2))
return g / g.sum()
def __init__(self, im, sigma_spatial=12, sigma_luma=4, sigma_chroma=4):
im_yuv = rgb2yuv(im)
# Compute 5-dimensional XYLUV bilateral-space coordinates
Iy, Ix = np.mgrid[:im.shape[0], :im.shape[1]]
x_coords = (Ix / sigma_spatial).astype(int)
y_coords = (Iy / sigma_spatial).astype(int)
luma_coords = (im_yuv[..., 0] /sigma_luma).astype(int)
chroma_coords = (im_yuv[..., 1:] / sigma_chroma).astype(int)
coords = np.dstack((x_coords, y_coords, luma_coords, chroma_coords))
coords_flat = coords.reshape(-1, coords.shape[-1])
self.npixels, self.dim = coords_flat.shape
# Hacky "hash vector" for coordinates,
# Requires all scaled coordinates be < MAX_VAL
self.hash_vec = (MAX_VAL**np.arange(self.dim))
# Construct S and B matrix
self._compute_factorization(coords_flat)
def get_local_mesh(self):
# Create the mesh
X = np.mgrid[self.rank*self.Np[0]:(self.rank+1)*self.Np[0], :self.N[1]].astype(self.float)
X[0] *= self.L[0]/self.N[0]
X[1] *= self.L[1]/self.N[1]
return X
def get_local_mesh(self):
xyrank = self.comm0.Get_rank() # Local rank in xz-plane
yzrank = self.comm1.Get_rank() # Local rank in xy-plane
# Create the physical mesh
x1 = slice(xyrank * self.N1[0], (xyrank+1) * self.N1[0], 1)
x2 = slice(yzrank * self.N2[1], (yzrank+1) * self.N2[1], 1)
X = np.mgrid[x1, x2, :self.N[2]].astype(self.float)
X[0] *= self.L[0]/self.N[0]
X[1] *= self.L[1]/self.N[1]
X[2] *= self.L[2]/self.N[2]
return X
def get_tform_coords(im_size):
coords0, coords1, coords2 = np.mgrid[:im_size[0], :im_size[1], :im_size[2]]
coords = np.array([coords0 - im_size[0] / 2, coords1 - im_size[1] / 2, coords2 - im_size[2] / 2])
return np.append(coords.reshape(3, -1), np.ones((1, np.prod(im_size))), axis=0)
def _FSpecialGauss(size, sigma):
"""Function to mimic the 'fspecial' gaussian MATLAB function."""
radius = size // 2
offset = 0.0
start, stop = -radius, radius + 1
if size % 2 == 0:
offset = 0.5
stop -= 1
x, y = np.mgrid[offset + start:stop, offset + start:stop]
assert len(x) == size
g = np.exp(-((x ** 2 + y ** 2) / (2.0 * sigma ** 2)))
return g / g.sum()
def replace_field(f, mask):
"""Interpolates positions in field according to mask with a 2D cubic interpolator"""
lx, ly = f.shape
x, y = np.mgrid[0:lx, 0:ly]
C = CT_intp((x[~mask],y[~mask]),f[~mask], fill_value=0)
return C(x, y)
def get_frame(self, i, j):
"""
Perform interpolation to produce the deformed window for correlation.
This function takes the previously set displacement and interpolates the image for these coordinates.
If the cubic interpolation method is chosen, the cubic interpolation of this API is use.
For the bilinear method the build in scipy method `map_coordinates <https://goo.gl/wucmUO>`_ is used with *order* set to 1.
:param int i: first index in grid coordinates
:param int j: second index in grid coordinates
:returns: interpolated window for the grid coordinates i,j and the image set in initialization
"""
dws = self._shape[-1]
offset_x, offset_y = np.mgrid[-dws/2+0.5:dws/2+0.5, -dws/2+0.5:dws/2+0.5]
gx, gy = np.mgrid[0:dws, 0:dws]
grid_x = gx + self._distance*i
grid_y = gy + self._distance*j
ptsax = (grid_x + self._u_disp(i, j, offset_x, offset_y)).ravel()
ptsay = (grid_y + self._v_disp(i, j, offset_x, offset_y)).ravel()
p, q = self._shape[-2:]
if self._ipmethod == 'bilinear':
return map_coordinates(self._frame, [ptsax, ptsay], order=1).reshape(p, q)
if self._ipmethod == 'cubic':
return self._cube_ip.interpolate(ptsax, ptsay).reshape(p, q)
def makeMTX(spat_coeffs, radial_filter, kr_IDX, viz_order=None, stepsize_deg=1):
"""Returns a plane wave decomposition over a full sphere
Parameters
----------
spat_coeffs : array_like
Spatial fourier coefficients
radial_filter : array_like
Modal radial filters
kr_IDX : int
Index of kr to be computed
viz_order : int, optional
Order of the spatial fourier transform [Default: Highest available]
stepsize_deg : float, optional
Integer Factor to increase the resolution. [Default: 1]
Returns
-------
mtxData : array_like
Plane wave decomposition (frequency domain)
Note
----
The file generates a Matrix of 181x360 pixels for the
visualisation with visualize3D() in 1[deg] Steps (65160 plane waves).
"""
if not viz_order:
viz_order = _np.int(_np.ceil(_np.sqrt(spat_coeffs.shape[0]) - 1))
angles = _np.mgrid[0:360:stepsize_deg, 0:181:stepsize_deg].reshape((2, -1)) * _np.pi / 180
Y = plane_wave_decomp(viz_order, angles, spat_coeffs[:, kr_IDX], radial_filter[:, kr_IDX])
return Y.reshape((360, -1)).T # Return pwd data as [181, 360] matrix
def plot3Dgrid(rows, cols, viz_data, style, normalize=True, title=None):
if len(viz_data) > rows * cols:
raise ValueError('Number of plot data is more than the specified rows and columns.')
fig = tools.make_subplots(rows, cols, specs=[[{'is_3d': True}] * cols] * rows, print_grid=False)
if style == 'flat':
layout_3D = dict(
xaxis=dict(range=[0, 360]),
yaxis=dict(range=[0, 181]),
aspectmode='manual',
aspectratio=dict(x=3.6, y=1.81, z=1)
)
else:
layout_3D = dict(
xaxis=dict(range=[-1, 1]),
yaxis=dict(range=[-1, 1]),
zaxis=dict(range=[-1, 1]),
aspectmode='cube'
)
rows, cols = _np.mgrid[1:rows + 1, 1: cols + 1]
rows = rows.flatten()
cols = cols.flatten()
for IDX in range(0, len(viz_data)):
cur_row = rows[IDX]
cur_col = cols[IDX]
fig.append_trace(genVisual(viz_data[IDX], style=style, normalize=normalize), cur_row, cur_col)
fig.layout['scene' + str(IDX + 1)].update(layout_3D)
if title is not None:
fig.layout.update(title=title)
filename = title + '.html'
else:
filename = str(current_time()) + '.html'
if env_info() == 'jupyter_notebook':
plotly_off.iplot(fig)
else:
plotly_off.plot(fig, filename=filename)
def gk(c1,r1,c2,r2):
# First, create X and Y arrays indicating distance to the boundaries of the paintbrush
# In this current context, im is the ordinal number of pixels (64 typically)
sigma = 0.3
im = 64
x = np.repeat([np.concatenate([np.mgrid[-c1:0],np.zeros(c2-c1),np.mgrid[1:1+im-c2]])],im,axis=0)
y = np.repeat(np.vstack(np.concatenate([np.mgrid[-r1:0],np.zeros(r2-r1),np.mgrid[1:1+im-r2]])),im,axis=1)
g = np.exp(-(x**2/float(im)+y**2/float(im))/(2*sigma**2))
return np.repeat([g],3,axis=0) # remove the 3 if you want to apply this to mask rather than an RGB channel
# This function reduces the likelihood of a change based on how close each individual pixel is to a maximal value.
# Consider conditioning this based on the gK value and the requested color. I.E. instead of just a flat distance from 128,
# have it be a difference from the expected color at a given location. This could also be used to "weight" the image towards staying the same.
def fcn_FDEM_InductionSpherePlaneWidget(xtx,ytx,ztx,m,orient,x0,y0,z0,a,sig,mur,xrx,yrx,zrx,logf,Comp,Phase):
sig = 10**sig
f = 10**logf
fvec = np.logspace(0,8,41)
xmin, xmax, dx, ymin, ymax, dy = -30., 30., 0.3, -30., 30., 0.4
X,Y = np.mgrid[xmin:xmax+dx:dx, ymin:ymax+dy:dy]
X = np.transpose(X)
Y = np.transpose(Y)
Obj = SphereFEM(m,orient,xtx,ytx,ztx)
Hx,Hy,Hz,Habs = Obj.fcn_ComputeFrequencyResponse(f,sig,mur,a,x0,y0,z0,X,Y,zrx)
Hxi,Hyi,Hzi,Habsi = Obj.fcn_ComputeFrequencyResponse(fvec,sig,mur,a,x0,y0,z0,xrx,yrx,zrx)
fig1 = plt.figure(figsize=(17,6))
Ax1 = fig1.add_axes([0.04,0,0.43,1])
Ax2 = fig1.add_axes([0.6,0,0.4,1])
if Comp == 'x':
Ax1 = plotAnomalyXYplane(Ax1,f,X,Y,ztx,Hx,Comp,Phase)
Ax1 = plotPlaceTxRxSphereXY(Ax1,xtx,ytx,xrx,yrx,x0,y0,a)
Ax2 = plotResponseFEM(Ax2,f,fvec,Hxi,Comp)
elif Comp == 'y':
Ax1 = plotAnomalyXYplane(Ax1,f,X,Y,ztx,Hy,Comp,Phase)
Ax1 = plotPlaceTxRxSphereXY(Ax1,xtx,ytx,xrx,yrx,x0,y0,a)
Ax2 = plotResponseFEM(Ax2,f,fvec,Hyi,Comp)
elif Comp == 'z':
Ax1 = plotAnomalyXYplane(Ax1,f,X,Y,ztx,Hz,Comp,Phase)
Ax1 = plotPlaceTxRxSphereXY(Ax1,xtx,ytx,xrx,yrx,x0,y0,a)
Ax2 = plotResponseFEM(Ax2,f,fvec,Hzi,Comp)
elif Comp == 'abs':
Ax1 = plotAnomalyXYplane(Ax1,f,X,Y,ztx,Habs,Comp,Phase)
Ax1 = plotPlaceTxRxSphereXY(Ax1,xtx,ytx,xrx,yrx,x0,y0,a)
Ax2 = plotResponseFEM(Ax2,f,fvec,Habsi,Comp)
plt.show(fig1)
def rebin(a, newshape):
"""Rebin an array to a new shape."""
assert len(a.shape) == len(newshape)
slices = [slice(0, old, float(old) / new)
for old, new in zip(a.shape, newshape)]
coordinates = np.mgrid[slices]
indices = coordinates.astype('i')
return a[tuple(indices)]
