def get_mag_ang(img):
"""
Gets image gradient (magnitude) and orientation (angle)
Args:
img
Returns:
Gradient, orientation
"""
img = np.sqrt(img)
gx = cv2.Sobel(np.float32(img), cv2.CV_32F, 1, 0)
gy = cv2.Sobel(np.float32(img), cv2.CV_32F, 0, 1)
mag, ang = cv2.cartToPolar(gx, gy)
return mag, ang, gx, gy
python类gradient()的实例源码
def calculate_angle_quality_thread(self,
voi,
gantry,
couch,
calculate_from=0,
stepsize=1.0,
q=None,
avoid=[],
voi_cube=None,
gradient=True):
""" TODO: Documentation
"""
os.nice(1)
for couch_angle in couch:
qual = self.calculate_angle_quality(voi, gantry, couch_angle, calculate_from, stepsize, avoid, voi_cube,
gradient)
q.put({"couch": couch_angle, "gantry": gantry, "data": qual})
def test_model_monovariate(func, var, par, k):
model = Model(func, var, par)
x, dx = np.linspace(0, 10, 100, retstep=True, endpoint=False)
U = np.cos(x * 2 * np.pi / 10)
fields = model.fields_template(x=x, U=U)
parameters = dict(periodic=True, k=k)
F = model.F(fields, parameters)
J_sparse = model.J(fields, parameters)
J_dense = model.J(fields, parameters, sparse=False)
J_approx = model.F.diff_approx(fields, parameters)
dxU = np.gradient(np.pad(U, 2, mode='wrap')) / dx
dxxU = np.gradient(dxU) / dx
dxU = dxU[2:-2]
dxxU = dxxU[2:-2]
assert np.isclose(F, k * dxxU, atol=1E-3).all()
assert np.isclose(J_approx, J_sparse.todense(), atol=1E-3).all()
assert np.isclose(J_approx, J_dense, atol=1E-3).all()
def test_model_bivariate(func, var, par):
model = Model(func, var, par)
x, dx = np.linspace(0, 10, 100, retstep=True, endpoint=False)
u = np.cos(x * 2 * np.pi / 10)
v = np.sin(x * 2 * np.pi / 10)
fields = model.fields_template(x=x, u=u, v=v)
parameters = dict(periodic=True, k1=1, k2=1)
F = model.F(fields, parameters)
J_sparse = model.J(fields, parameters)
J_dense = model.J(fields, parameters, sparse=False)
J_approx = model.F.diff_approx(fields, parameters)
dxu = np.gradient(np.pad(u, 2, mode='wrap')) / dx
dxxu = np.gradient(dxu) / dx
dxu = dxu[2:-2]
dxxu = dxxu[2:-2]
dxv = np.gradient(np.pad(v, 2, mode='wrap')) / dx
dxxv = np.gradient(dxv) / dx
dxv = dxv[2:-2]
dxxv = dxxv[2:-2]
assert np.isclose(F, np.vstack([dxxv, dxxu]).flatten('F'),
atol=1E-3).all()
assert np.isclose(J_approx, J_sparse.todense(), atol=1E-3).all()
assert np.isclose(J_approx, J_dense, atol=1E-3).all()
def inflection_points(points, axis, span):
'''
Find the list of vertices that preceed inflection points in a curve. The curve is differentiated
with respect to the coordinate system defined by axis and span.
axis: A vector representing the vertical axis of the coordinate system.
span: A vector representing the the horiztonal axis of the coordinate system.
returns: a list of points in space corresponding to the vertices that
immediately preceed inflection points in the curve
'''
coords_on_span = points.dot(span)
dx = np.gradient(coords_on_span)
coords_on_axis = points.dot(axis)
# Take the second order finite difference of the curve with respect to the
# defined coordinate system
finite_difference_2 = np.gradient(np.gradient(coords_on_axis, dx), dx)
# Compare the product of all neighboring pairs of points in the second derivative
# If a pair of points has a negative product, then the second derivative changes sign
# at one of those points, signalling an inflection point
is_inflection_point = [finite_difference_2[i] * finite_difference_2[i + 1] <= 0 for i in range(len(finite_difference_2) - 1)]
inflection_point_indices = [i for i, b in enumerate(is_inflection_point) if b]
if len(inflection_point_indices) == 0: # pylint: disable=len-as-condition
return []
return points[inflection_point_indices]
def calculate_sift_grid(self, image, gridH, gridW):
H, W = image.shape
Npatches = gridH.size
feaArr = np.zeros((Npatches, Nsamples * Nangles))
# calculate gradient
GH, GW = gen_dgauss(self.sigma)
IH = signal.convolve2d(image, GH, mode='same')
IW = signal.convolve2d(image, GW, mode='same')
Imag = np.sqrt(IH ** 2 + IW ** 2)
Itheta = np.arctan2(IH, IW)
Iorient = np.zeros((Nangles, H, W))
for i in range(Nangles):
Iorient[i] = Imag * np.maximum(np.cos(Itheta - angles[i]) ** alpha, 0)
for i in range(Npatches):
currFeature = np.zeros((Nangles, Nsamples))
for j in range(Nangles):
currFeature[j] = np.dot(self.weights, \
Iorient[j, gridH[i]:gridH[i] + self.pS, gridW[i]:gridW[i] + self.pS].flatten())
feaArr[i] = currFeature.flatten()
return feaArr
def extract_sift_patches(self, image, gridH, gridW):
# extracts the sift descriptor of patches
# in positions (gridH,gridW) in the image
H, W = image.shape
Npatches = gridH.size
feaArr = np.zeros((Npatches, Nsamples * Nangles))
# calculate gradient
GH, GW = gen_dgauss(self.sigma)
IH = signal.convolve2d(image, GH, mode='same')
IW = signal.convolve2d(image, GW, mode='same')
Imag = np.sqrt(IH ** 2 + IW ** 2)
Itheta = np.arctan2(IH, IW)
Iorient = np.zeros((Nangles, H, W))
for i in range(Nangles):
Iorient[i] = Imag * np.maximum(np.cos(Itheta - angles[i]) ** alpha, 0)
for i in range(Npatches):
currFeature = np.zeros((Nangles, Nsamples))
for j in range(Nangles):
currFeature[j] = np.dot(self.weights, \
Iorient[j, gridH[i]:gridH[i] + self.pS, gridW[i]:gridW[i] + self.pS].flatten())
feaArr[i] = currFeature.flatten()
# feaArr contains each descriptor in a row
feaArr = self.normalize_sift(feaArr)
return feaArr
def test_specific_axes(self):
# Testing that gradient can work on a given axis only
v = [[1, 1], [3, 4]]
x = np.array(v)
dx = [np.array([[2., 3.], [2., 3.]]),
np.array([[0., 0.], [1., 1.]])]
assert_array_equal(gradient(x, axis=0), dx[0])
assert_array_equal(gradient(x, axis=1), dx[1])
assert_array_equal(gradient(x, axis=-1), dx[1])
assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])
# test axis=None which means all axes
assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
# and is the same as no axis keyword given
assert_almost_equal(gradient(x, axis=None), gradient(x))
# test vararg order
assert_array_equal(gradient(x, 2, 3, axis=(1, 0)), [dx[1]/2.0, dx[0]/3.0])
# test maximal number of varargs
assert_raises(SyntaxError, gradient, x, 1, 2, axis=1)
assert_raises(ValueError, gradient, x, axis=3)
assert_raises(ValueError, gradient, x, axis=-3)
assert_raises(TypeError, gradient, x, axis=[1,])
def contour_from_array_at_label_l(im_arr, l, thr=0.3, omit_axis=None, verbose=0):
if verbose > 0:
print(l)
array_label_l = im_arr == l
assert isinstance(array_label_l, np.ndarray)
gra = np.gradient(array_label_l.astype(np.bool).astype(np.float64))
if omit_axis is None:
norm_gra = np.sqrt(gra[0] ** 2 + gra[1] ** 2 + gra[2] ** 2) > thr
elif omit_axis == 'x':
norm_gra = np.sqrt(gra[1] ** 2 + gra[2] ** 2) > thr
elif omit_axis == 'y':
norm_gra = np.sqrt(gra[0] ** 2 + gra[2] ** 2) > thr
elif omit_axis == 'z':
norm_gra = np.sqrt(gra[0] ** 2 + gra[1] ** 2) > thr
else:
raise IOError
return norm_gra
def computeWeightsLocallyNormalized(I, centered_gradient=True, norm_radius=45):
h,w = I.shape[:2]
if centered_gradient:
gy,gx = np.gradient(I)[:2]
gysq = (gy**2).mean(axis=2) if gy.ndim > 2 else gy**2
gxsq = (gx**2).mean(axis=2) if gx.ndim > 2 else gx**2
gxsq_local_mean = cv2.blur(gxsq, ksize=(norm_radius, norm_radius))
gysq_local_mean = cv2.blur(gysq, ksize=(norm_radius, norm_radius))
w_horizontal = np.exp( - gxsq * 1.0/(2*np.maximum(1e-6, gxsq_local_mean)))
w_vertical = np.exp( - gysq * 1.0/(2*np.maximum(1e-6, gysq_local_mean)))
else:
raise Exception("NotImplementedYet")
return w_horizontal, w_vertical
def rel_vort(u,v,dx,dy):
#get the gradient of the wind in the u- and v-directions
du = np.gradient(u)
dv = np.gradient(v)
#compute the relative vorticity (units : 10^-5 s^-1)
vort = ((dv[-1]/dx) - (du[-2]/dy))*pow(10,5)
#return the vorticity (Units: 10^-5 s-1)
return vort
###############################################################################
########################## End function rel_vort() ##########################
###############################################################################
#################### 25.Begin function of wind_chill() ######################
## Required libraries: numpy #
## #
## Inputs: t2m = ndarray of 2-meter temperature values (Units: F) #
## #
## wind = ndarray of 10-meter bulk wind values (Units: MPH) #
## #
###############################################################################
def richardson(self):
u = self.dataSet.readNCVariable('U')
v = self.dataSet.readNCVariable('V')
u_corr = wrf.unstaggerX(u)
v_corr = wrf.unstaggerY(v)
wind = (u_corr*u_corr + v_corr*v_corr)**(0.5)
height = wrf.unstaggerZ(self.height)
#compute the vertical gradient of theta
dtheta = np.gradient(self.theta)[0]
du = np.gradient(wind)[0]
dz = np.gradient(height)[0]
#compute the richardson number
self.u10 = u_corr
self.v10 = v_corr
self.var = ((self.g/self.theta)*(dtheta/dz))/(du/dz)**(2)
self.var2 = self.var
self.varTitle = "Richardson Number\n" + self.dataSet.getTime()
#Set short variable title for time series
self.sTitle = "Richardson Number"
def smooth_abs(x, dx=0.01):
"""smoothed version of absolute vaue function, with quadratic instead of sharp bottom.
Derivative w.r.t. variable of interest. Width of quadratic can be controlled"""
x, n = _checkIfFloat(x)
y = np.abs(x)
idx = np.logical_and(x > -dx, x < dx)
y[idx] = x[idx]**2/(2.0*dx) + dx/2.0
# gradient
dydx = np.ones_like(x)
dydx[x <= -dx] = -1.0
dydx[idx] = x[idx]/dx
if n == 1:
y = y[0]
dydx = dydx[0]
return y, dydx
def test_specific_axes(self):
# Testing that gradient can work on a given axis only
v = [[1, 1], [3, 4]]
x = np.array(v)
dx = [np.array([[2., 3.], [2., 3.]]),
np.array([[0., 0.], [1., 1.]])]
assert_array_equal(gradient(x, axis=0), dx[0])
assert_array_equal(gradient(x, axis=1), dx[1])
assert_array_equal(gradient(x, axis=-1), dx[1])
assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])
# test axis=None which means all axes
assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
# and is the same as no axis keyword given
assert_almost_equal(gradient(x, axis=None), gradient(x))
# test vararg order
assert_array_equal(gradient(x, 2, 3, axis=(1, 0)), [dx[1]/2.0, dx[0]/3.0])
# test maximal number of varargs
assert_raises(SyntaxError, gradient, x, 1, 2, axis=1)
assert_raises(ValueError, gradient, x, axis=3)
assert_raises(ValueError, gradient, x, axis=-3)
assert_raises(TypeError, gradient, x, axis=[1,])
def main():
"""Load image, compute derivatives, plot."""
img = data.camera()
imgy_np, imgx_np = np.gradient(img)
imgx_ski, imgy_ski = _compute_derivatives(img)
# dx
util.plot_images_grayscale(
[img, dx_symmetric(img), dx_forward(img), dx_backward(img), imgx_np, imgx_ski],
["Image", "dx (symmetric)", "dx (forward)", "dx (backward)",
"Ground Truth (numpy)", "Ground Truth (scikit-image)"]
)
# dy
util.plot_images_grayscale(
[img, dy_symmetric(img), dy_forward(img), dy_backward(img), imgy_np, imgy_ski],
["Image", "dy (symmetric)", "dy (forward)", "dy (backward)",
"Ground Truth (numpy)", "Ground Truth (scikit-image)"]
)
def test_specific_axes(self):
# Testing that gradient can work on a given axis only
v = [[1, 1], [3, 4]]
x = np.array(v)
dx = [np.array([[2., 3.], [2., 3.]]),
np.array([[0., 0.], [1., 1.]])]
assert_array_equal(gradient(x, axis=0), dx[0])
assert_array_equal(gradient(x, axis=1), dx[1])
assert_array_equal(gradient(x, axis=-1), dx[1])
assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])
# test axis=None which means all axes
assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
# and is the same as no axis keyword given
assert_almost_equal(gradient(x, axis=None), gradient(x))
# test vararg order
assert_array_equal(gradient(x, 2, 3, axis=(1, 0)), [dx[1]/2.0, dx[0]/3.0])
# test maximal number of varargs
assert_raises(SyntaxError, gradient, x, 1, 2, axis=1)
assert_raises(ValueError, gradient, x, axis=3)
assert_raises(ValueError, gradient, x, axis=-3)
assert_raises(TypeError, gradient, x, axis=[1,])
def direction_map(dmap: DistanceMap):
r"""Computes normalized gradient of distance map. Not defined when length of
the gradient is zero.
.. math::
\hat{\mathbf{e}}_{S} = -\frac{\nabla S(\mathbf{x})}{\| \nabla S(\mathbf{x}) \|}
Args:
dmap (numpy.ndarray):
Distance map.
Returns:
DirectionMap: Direction map.
"""
u, v = np.gradient(dmap)
l = np.hypot(u, v)
# Avoids zero division
l[l == 0] = np.nan
# Flip order from (row, col) to (x, y)
return v / l, u / l
# Potentials
def __entropy(self, spectrum, m):
"""Calculates get m-th derivative of the real part of spectrum and
returns entropy of its absolute value. """
assert type(m) is int, 'm should be a (possitive) integer'
assert m > 0, 'need m > 0'
#get the real part of the spectrum
spect = np.array(spectrum)
spect = spect.real
# calculate the m-th derivative of the real part of the spectrum
spectrumDerivative = spect
for i in range(m):
spectrumDerivative = np.gradient(spectrumDerivative)
# now get the entropy of the abslolute value of the m-th derivative:
entropy = sp.stats.entropy(np.abs(spectrumDerivative))
return entropy
def threshold_gradients(self, grad_thresh):
"""Creates a new DepthImage by zeroing out all depths
where the magnitude of the gradient at that point is
greater than grad_thresh.
Parameters
----------
grad_thresh : float
A threshold for the gradient magnitude.
Returns
-------
:obj:`DepthImage`
A new DepthImage created from the thresholding operation.
"""
data = np.copy(self._data)
gx, gy = self.gradients()
gradients = np.zeros([gx.shape[0], gx.shape[1], 2])
gradients[:, :, 0] = gx
gradients[:, :, 1] = gy
gradient_mags = np.linalg.norm(gradients, axis=2)
ind = np.where(gradient_mags > grad_thresh)
data[ind[0], ind[1]] = 0.0
return DepthImage(data, self._frame)
def normal_cloud_im(self):
"""Generate a NormalCloudImage from the PointCloudImage.
Returns
-------
:obj:`NormalCloudImage`
The corresponding NormalCloudImage.
"""
gx, gy, _ = np.gradient(self.data)
gx_data = gx.reshape(self.height * self.width, 3)
gy_data = gy.reshape(self.height * self.width, 3)
pc_grads = np.cross(gx_data, gy_data) # default to point toward camera
pc_grad_norms = np.linalg.norm(pc_grads, axis=1)
normal_data = pc_grads / np.tile(pc_grad_norms[:, np.newaxis], [1, 3])
normal_im_data = normal_data.reshape(self.height, self.width, 3)
return NormalCloudImage(normal_im_data, frame=self.frame)
def visualize_derivatives(image):
'''
Plot gradient on left and Laplacian on right.
Only tested on 2D 1-channel float imags
'''
dx1,dy1 = np.gradient(image)
gradient = dx1 + 1j*dy1
a1 = np.abs(gradient)
plt.figure(None,(12,6))
plt.subplot(121)
a1 = mean_center(blur(exposure.equalize_hist(unitize(a1)),1))
plt.imshow(a1,
origin='lower',interpolation='nearest',cmap='gray',extent=(0,64,)*2)
plt.title('Gradient Magnitude')
plt.subplot(122)
laplacian = scipy.ndimage.filters.laplace(image)
lhist = mean_center(
blur(exposure.equalize_hist(unitize(laplacian)),1))
plt.imshow(lhist,
origin='lower',interpolation='nearest',cmap='gray',extent=(0,64,)*2)
plt.title('Laplacian')
return gradient, laplacian
def test_specific_axes(self):
# Testing that gradient can work on a given axis only
v = [[1, 1], [3, 4]]
x = np.array(v)
dx = [np.array([[2., 3.], [2., 3.]]),
np.array([[0., 0.], [1., 1.]])]
assert_array_equal(gradient(x, axis=0), dx[0])
assert_array_equal(gradient(x, axis=1), dx[1])
assert_array_equal(gradient(x, axis=-1), dx[1])
assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])
# test axis=None which means all axes
assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
# and is the same as no axis keyword given
assert_almost_equal(gradient(x, axis=None), gradient(x))
# test vararg order
assert_array_equal(gradient(x, 2, 3, axis=(1, 0)), [dx[1]/2.0, dx[0]/3.0])
# test maximal number of varargs
assert_raises(SyntaxError, gradient, x, 1, 2, axis=1)
assert_raises(ValueError, gradient, x, axis=3)
assert_raises(ValueError, gradient, x, axis=-3)
assert_raises(TypeError, gradient, x, axis=[1,])
def _get_displacement_function(self, f):
"""
Getter function for calculating the displacement.
:param f: field that is used for the displacement, mainly velocity components
:returns: function of the Taylor expanded field to first order
"""
dx = self._distance
f_x, f_y = np.gradient(f , dx)
f_xx, f_xy = np.gradient(f_x, dx)
f_yx, f_yy = np.gradient(f_y, dx)
return lambda i, j, x, y : (f[i, j] + x*f_x[i, j] + y*f_y[i, j]
+ 0.5*(f_xx[i, j]*x**2 + 2*f_xy[i, j]*x*y + f_yy[i, j]*y**2))
#For the bilinear method the build in scipy method `map_coordinates <https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.map_coordinates.html>`_ is used with *order* set to 1.
def psf_shift(lenslet_efl, dx, dy, mag=1):
''' Computes the shift of a PSF, in microns.
Args:
lenslet_efl (`microns`): EFL of lenslets.
dx (`np.ndarray`): dx gradient of wavefront.
dy (`np.ndarray`): dy gradient of wavefront.
mag (`float`): magnification of the collimation system.
Returns:
`numpy.ndarray`: array of PSF shifts.
Notes:
see eq. 12 of Chanan, "Principles of Wavefront Sensing and Reconstruction"
delta = m * fl * grad(z)
m is magnification of SH system
fl is lenslet focal length
grad(z) is the x, y gradient of the opd, z, which is expressed in
radians.
'''
coef = -mag * lenslet_efl
return coef * dx, coef * dy
def test_basic(self):
v = [[1, 1], [3, 4]]
x = np.array(v)
dx = [np.array([[2., 3.], [2., 3.]]),
np.array([[0., 0.], [1., 1.]])]
assert_array_equal(gradient(x), dx)
assert_array_equal(gradient(v), dx)
def test_badargs(self):
# for 2D array, gradient can take 0, 1, or 2 extra args
x = np.array([[1, 1], [3, 4]])
assert_raises(SyntaxError, gradient, x, np.array([1., 1.]),
np.array([1., 1.]), np.array([1., 1.]))
def test_masked(self):
# Make sure that gradient supports subclasses like masked arrays
x = np.ma.array([[1, 1], [3, 4]],
mask=[[False, False], [False, False]])
out = gradient(x)[0]
assert_equal(type(out), type(x))
# And make sure that the output and input don't have aliased mask
# arrays
assert_(x.mask is not out.mask)
# Also check that edge_order=2 doesn't alter the original mask
x2 = np.ma.arange(5)
x2[2] = np.ma.masked
np.gradient(x2, edge_order=2)
assert_array_equal(x2.mask, [False, False, True, False, False])
def test_datetime64(self):
# Make sure gradient() can handle special types like datetime64
x = np.array(
['1910-08-16', '1910-08-11', '1910-08-10', '1910-08-12',
'1910-10-12', '1910-12-12', '1912-12-12'],
dtype='datetime64[D]')
dx = np.array(
[-5, -3, 0, 31, 61, 396, 731],
dtype='timedelta64[D]')
assert_array_equal(gradient(x), dx)
assert_(dx.dtype == np.dtype('timedelta64[D]'))
def test_timedelta64(self):
# Make sure gradient() can handle special types like timedelta64
x = np.array(
[-5, -3, 10, 12, 61, 321, 300],
dtype='timedelta64[D]')
dx = np.array(
[2, 7, 7, 25, 154, 119, -21],
dtype='timedelta64[D]')
assert_array_equal(gradient(x), dx)
assert_(dx.dtype == np.dtype('timedelta64[D]'))
def test_specific_axes(self):
# Testing that gradient can work on a given axis only
v = [[1, 1], [3, 4]]
x = np.array(v)
dx = [np.array([[2., 3.], [2., 3.]]),
np.array([[0., 0.], [1., 1.]])]
assert_array_equal(gradient(x, axis=0), dx[0])
assert_array_equal(gradient(x, axis=1), dx[1])
assert_array_equal(gradient(x, axis=-1), dx[1])
assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]])
# test axis=None which means all axes
assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]])
# and is the same as no axis keyword given
assert_almost_equal(gradient(x, axis=None), gradient(x))
# test vararg order
assert_array_equal(gradient(x, 2, 3, axis=(1, 0)), [dx[1]/2.0, dx[0]/3.0])
# test maximal number of varargs
assert_raises(SyntaxError, gradient, x, 1, 2, axis=1)
assert_raises(ValueError, gradient, x, axis=3)
assert_raises(ValueError, gradient, x, axis=-3)
assert_raises(TypeError, gradient, x, axis=[1,])
