Python Numpy gradient原始碼解析
阿新 • • 發佈：2018-12-21
複習影象梯度，發現Numpy有一個梯度計算函式，解析它的原始碼和需要注意的問題，最後自定義一個梯度函式
目錄
用法解析
Numpy提供了陣列梯度計算函式
gradient(f, *varargs, **kwargs)
輸入
- 必選引數：類N維陣列（列表/元組/陣列）
- 可選引數：標量列表或陣列列表，用於計算差分時的間隔空間
- 單個標量：為所有軸指定間隔
N個標量：為每一軸指定了間隔，比如dx, dx, dz...N個數組：為f的每一軸都指定了取樣值的下標，每個陣列的長度必須匹配相對應那一維的長度N個標量或陣列的組合，其作用等同於2和3- 注意：如果關鍵字引數
axis被指定，那麼這個可選引數的長度必須和要計算的軸數一致
- 關鍵字引數有
2個：edge_order：在邊界使用第N-階差分公式，預設是1-階axis：None或者是整數又或者是整數的元組，表示沿著指定的軸計算梯度，預設為None，表示計算所有軸的梯度。如果輸入為負數，表示從後向前計算軸的梯度
從它的介紹中，我唯一沒理解的就是可選引數的用法，看了程式碼後也沒理解，因為它和一階差分沒有關係，所以在後面的內容就忽略對可選引數的解析
輸出
Return the gradient of an N-dimensional array. The returned gradient hence has the same shape as the input array.
如果輸入是一維陣列，返回同樣大小的梯度陣列
如果是多維陣列，並要求計算多維，那麼返回一個列表，包含計算的每一維的梯度，其大小和輸入陣列一致
計算方式
The gradient is computed using second order accurate central differences
in the interior points and either first or second order accurate one-sides
(forward or backwards) differences at the boundaries.
在它的介紹中，中間元素使用**2階中心差分**，但是檢視原始碼發現，它預設使用的其實是一階中心差分
3和4選項，但是程式碼中的二階差分公式也沒看懂？，應該是使用泰勒級數的近似公式，所以跳過對程式碼中二階差分的介紹
兩邊使用一階前向或後向差分，也可通過關鍵字引數edge_order修改計算方式
差分公式
一階前向差分：
一階後向差分：
一階中心差分：
二階中心差分：
示例和問題
示例和問題
計算一維或二維陣列的梯度，以及計算影象梯度
一維陣列
>>> import numpy as np
>>> arr1 = np.array(range(5))
>>> arr1
array([0, 1, 2, 3, 4])
>>> np.gradient(arr1)
array([1., 1., 1., 1., 1.])
>>> np.gradient(arr1, 2) # 增加可選引數
array([0.5, 0.5, 0.5, 0.5, 0.5])
問題：沒理解可選引數的適用場景？
二維陣列
>>> arr2 = np.random.randint(1, 9, (3,4))
>>> arr2
array([[6, 3, 5, 1],
[3, 8, 6, 3],
[7, 3, 3, 4]])
>>> np.gradient(arr2)
[array([[-3. , 5. , 1. , 2. ],
[ 0.5, 0. , -1. , 1.5],
[ 4. , -5. , -3. , 1. ]]), array([[-3. , -0.5, -1. , -4. ],
[ 5. , 1.5, -2.5, -3. ],
[-4. , -2. , 0.5, 1. ]])]
>>> np.gradient(arr2, axis=0) # 對行求導
array([[-3. , 5. , 1. , 2. ],
[ 0.5, 0. , -1. , 1.5],
[ 4. , -5. , -3. , 1. ]])
>>> np.gradient(arr2, axis=1) # 對列求導
array([[-3. , -0.5, -1. , -4. ],
[ 5. , 1.5, -2.5, -3. ],
[-4. , -2. , 0.5, 1. ]])
axis=0表示對行求導，axis=1表示對列求導
影象梯度計算需要注意的問題
>>> import cv2 as cv
>>> img = cv.imread("C:\\python\\ConvolveAndGradient\\lena.jpg", 0)
>>> arr3 = img[:3, :4]
>>> arr3
array([[163, 162, 161, 160],
[162, 162, 162, 161],
[162, 162, 163, 161]], dtype=uint8)
>>> np.gradient(arr3, axis=0)
array([[255. , 0. , 1. , 1. ],
[127.5, 0. , 1. , 0.5],
[ 0. , 0. , 1. , 0. ]])
>>> np.gradient(arr3, axis=1)
array([[255. , 127. , 127. , 255. ],
[ 0. , 0. , 127.5, 255. ],
[ 0. , 0.5, 127.5, 254. ]])
輸入影象計算梯度時，發現數組邊界的梯度計算有誤，後來看原始碼的時候突然發覺是因為數值型別的關係，影象型別是uint8，得到負數會溢位，所以需要先轉換成更高精度的型別
>>> arr4 = np.asarray(img[:3, :4], dtype=np.double)
>>> arr4
array([[163., 162., 161., 160.],
[162., 162., 162., 161.],
[162., 162., 163., 161.]])
>>> np.gradient(arr4, axis=0)
array([[-1. , 0. , 1. , 1. ],
[-0.5, 0. , 1. , 0.5],
[ 0. , 0. , 1. , 0. ]])
>>> np.gradient(arr4, axis=1)
array([[-1. , -1. , -1. , -1. ],
[ 0. , 0. , -0.5, -1. ],
[ 0. , 0.5, -0.5, -2. ]])
原始碼解析
def gradient(f, *varargs, **kwargs):
# 將類陣列轉換成陣列
f = np.asanyarray(f)
# 獲取維數
N = f.ndim # number of dimensions
# 有沒有指定計算的軸數，預設為空
axes = kwargs.pop('axis', None)
if axes is None:
# 如果為空，那麼計算每一軸的梯度
axes = tuple(range(N))
else:
# 如果不為空，將其轉換成非負整數軸數列表，比如，axis=(-2,-1) -> (0, 1)
axes = _nx.normalize_axis_tuple(axes, N)
# 待計算的軸數
len_axes = len(axes)
# 可選引數長度，不解析，預設為0
n = len(varargs)
if n == 0:
# no spacing argument - use 1 in all axes
dx = [1.0] * len_axes
elif n == 1 and np.ndim(varargs[0]) == 0:
# single scalar for all axes
dx = varargs * len_axes
elif n == len_axes:
# scalar or 1d array for each axis
dx = list(varargs)
for i, distances in enumerate(dx):
if np.ndim(distances) == 0:
continue
elif np.ndim(distances) != 1:
raise ValueError("distances must be either scalars or 1d")
if len(distances) != f.shape[axes[i]]:
raise ValueError("when 1d, distances must match "
"the length of the corresponding dimension")
diffx = np.diff(distances)
# if distances are constant reduce to the scalar case
# since it brings a consistent speedup
if (diffx == diffx[0]).all():
diffx = diffx[0]
dx[i] = diffx
else:
raise TypeError("invalid number of arguments")
# 邊界求導方式，預設一階，最多2階
edge_order = kwargs.pop('edge_order', 1)
if kwargs:
raise TypeError('"{}" are not valid keyword arguments.'.format(
'", "'.join(kwargs.keys())))
if edge_order > 2:
raise ValueError("'edge_order' greater than 2 not supported")
# use central differences on interior and one-sided differences on the
# endpoints. This preserves second order-accuracy over the full domain.
outvals = []
# create slice objects --- initially all are [:, :, ..., :]
# 為每一軸設定切片函式
slice1 = [slice(None)]*N
slice2 = [slice(None)]*N
slice3 = [slice(None)]*N
slice4 = [slice(None)]*N
otype = f.dtype
# 定義輸出數值型別，先解析輸入陣列數值型別
# 除特殊型別（np.datetime64/np.timedelta64/np.inexact）外，均轉換成np.double
if otype.type is np.datetime64:
# the timedelta dtype with the same unit information
otype = np.dtype(otype.name.replace('datetime', 'timedelta'))
# view as timedelta to allow addition
f = f.view(otype)
elif otype.type is np.timedelta64:
pass
elif np.issubdtype(otype, np.inexact):
pass
else:
# all other types convert to floating point
otype = np.double
# 求導每一軸
# 按axes定義的軸順序來求導
# 為啥要弄個dx？，它和可選引數有關，不解析可選引數，所以預設為1.0，不影響計算
for axis, ax_dx in zip(axes, dx):
if f.shape[axis] < edge_order + 1:
# 每一軸的長度必須符合邊界求導規則，比如一階求導，必須有2個數值；2階求導，必須有3個數值
raise ValueError(
"Shape of array too small to calculate a numerical gradient, "
"at least (edge_order + 1) elements are required.")
# result allocation
# 建立未初始化陣列
out = np.empty_like(f, dtype=otype)
# spacing for the current axis
# 因為不接戲可選引數，所以為True
uniform_spacing = np.ndim(ax_dx) == 0
# Numerical differentiation: 2nd order interior
# 切片1：取中間值
# 切片2：取前n-2個
# 切片3：取中間值
# 切片4：取後n-2個
slice1[axis] = slice(1, -1)
slice2[axis] = slice(None, -2)
slice3[axis] = slice(1, -1)
slice4[axis] = slice(2, None)
if uniform_spacing:
# 中間數值使用中心差分，有計算可知，是一階差分公式
out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / (2. * ax_dx)
else:
dx1 = ax_dx[0:-1]
dx2 = ax_dx[1:]
a = -(dx2)/(dx1 * (dx1 + dx2))
b = (dx2 - dx1) / (dx1 * dx2)
c = dx1 / (dx2 * (dx1 + dx2))
# fix the shape for broadcasting
shape = np.ones(N, dtype=int)
shape[axis] = -1
a.shape = b.shape = c.shape = shape
# 1D equivalent -- out[1:-1] = a * f[:-2] + b * f[1:-1] + c * f[2:]
out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
# Numerical differentiation: 1st order edges
if edge_order == 1:
# 前向差分
slice1[axis] = 0
slice2[axis] = 1
slice3[axis] = 0
dx_0 = ax_dx if uniform_spacing else ax_dx[0]
# 1D equivalent -- out[0] = (f[1] - f[0]) / (x[1] - x[0])
out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0
# 後向差分
slice1[axis] = -1
slice2[axis] = -1
slice3[axis] = -2
dx_n = ax_dx if uniform_spacing else ax_dx[-1]
# 1D equivalent -- out[-1] = (f[-1] - f[-2]) / (x[-1] - x[-2])
out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n
# Numerical differentiation: 2nd order edges
else:
slice1[axis] = 0
slice2[axis] = 0
slice3[axis] = 1
slice4[axis] = 2
if uniform_spacing:
a = -1.5 / ax_dx
b = 2. / ax_dx
c = -0.5 / ax_dx
else:
dx1 = ax_dx[0]
dx2 = ax_dx[1]
a = -(2. * dx1 + dx2)/(dx1 * (dx1 + dx2))
b = (dx1 + dx2) / (dx1 * dx2)
c = - dx1 / (dx2 * (dx1 + dx2))
# 1D equivalent -- out[0] = a * f[0] + b * f[1] + c * f[2]
out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
slice1[axis] = -1
slice2[axis] = -3
slice3[axis] = -2
slice4[axis] = -1
if uniform_spacing:
a = 0.5 / ax_dx
b = -2. / ax_dx
c = 1.5 / ax_dx
else:
dx1 = ax_dx[-2]
dx2 = ax_dx[-1]
a = (dx2) / (dx1 * (dx1 + dx2))
b = - (dx2 + dx1) / (dx1 * dx2)
c = (2. * dx2 + dx1) / (dx2 * (dx1 + dx2))
# 1D equivalent -- out[-1] = a * f[-3] + b * f[-2] + c * f[-1]
out[tuple(slice1)] = a * f[tuple(slice2)] + b * f[tuple(slice3)] + c * f[tuple(slice4)]
# 在結果陣列中新增該軸計算得到的梯度陣列
outvals.append(out)
# reset the slice object in this dimension to ":"
slice1[axis] = slice(None)
slice2[axis] = slice(None)
slice3[axis] = slice(None)
slice4[axis] = slice(None)
# 返回
if len_axes == 1:
return outvals[0]
else:
return outvals
自定義
下面自定義一個影象梯度計算函式
要求：
- 對輸入影象數值型別進行轉換，避免溢位
- 對中間畫素進行
1階中心差分，對邊界進行1階前向或後向差分 - 計算過程不改變輸入影象畫素值
- 返回一個列表，包含每一軸的梯度陣列
- 可以指定要計算哪一軸的梯度，如果僅計算一軸梯度，返回該陣列即可
實現如下：
def image_gradient(f, **kwargs):
"""
影象梯度求導
:param f: 類陣列
:return: 陣列或陣列列表，dtype=np.double
可選引數：axis - 空或整數或整數列表
"""
f = np.asanyarray(f, dtype=np.double)
N = f.ndim # number of dimensions
axes = kwargs.pop('axis', None)
if axes is None:
axes = tuple(range(N))
else:
axes = _nx.normalize_axis_tuple(axes, N)
len_axes = len(axes)
outvals = []
# create slice objects --- initially all are [:, :, ..., :]
slice1 = [slice(None)] * N
slice2 = [slice(None)] * N
slice3 = [slice(None)] * N
slice4 = [slice(None)] * N
otype = f.dtype
for axis in axes:
if f.shape[axis] < 2:
raise ValueError(
"Shape of array too small to calculate a numerical gradient, "
"at least 2 elements are required.")
# result allocation
out = np.empty_like(f, dtype=otype)
# Numerical differentiation: 2nd order interior
slice1[axis] = slice(1, -1)
slice2[axis] = slice(None, -2)
slice3[axis] = slice(1, -1)
slice4[axis] = slice(2, None)
# 1D equivalent -- out[0] = (f[1] - f[-1]) / 2
out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / 2
slice1[axis] = 0
slice2[axis] = 1
slice3[axis] = 0
# 1D equivalent -- out[0] = (f[1] - f[0]) / (x[1] - x[0])
out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)])
slice1[axis] = -1
slice2[axis] = -1
slice3[axis] = -2
# 1D equivalent -- out[-1] = (f[-1] - f[-2]) / (x[-1] - x[-2])
out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)])
outvals.append(out)
# reset the slice object in this dimension to ":"
slice1[axis] = slice(None)
slice2[axis] = slice(None)
slice3[axis] = slice(None)
slice4[axis] = slice(None)
if len_axes == 1:
return outvals[0]
else:
return outvals
輸出陣列或陣列列表的數值型別是np.double，計算影象梯度如下：
import cv2 as cv
import numpy as np
import numpy.core.numeric as _nx
import matplotlib.pyplot as plt
__author__ = 'zj'
image_name = 'lena.jpg'
def get_image():
return cv.imread(image_name, 0)
def show_image(grad_x, grad_y, grad):
plt.figure(figsize=(10, 5)) # 設定視窗大小
plt.suptitle('image_gradient') # 圖片名稱
plt.subplot(1, 3, 1)
plt.title('grad_x')
plt.imshow(grad_x, cmap='gray'), plt.axis('off')
plt.subplot(1, 3, 2)
plt.title('gray_y')
plt.imshow(grad_y, cmap='gray'), plt.axis('off')
plt.subplot(1, 3, 3)
plt.title('grad')
plt.imshow(grad, cmap='gray'), plt.axis('off')
plt.savefig('./grad.png') # 儲存影象
plt.show()
if __name__ == '__main__':
img = get_image()
grad_x = image_gradient(img, axis=1)
grad_y = image_gradient(img, axis=0)
grad = np.sqrt(grad_x ** 2 + grad_y ** 2)
abs_x = np.array(np.abs(grad_x), dtype=np.uint8)
abs_y = np.array(np.abs(grad_y), dtype=np.uint8)
abs_grad = np.array(np.abs(grad), dtype=np.uint8)
# cv.imshow("x", abs_x)
# cv.imshow("y", abs_y)
# cv.imshow("grad", abs_grad)
# cv.waitKey()
show_image(abs_x, abs_y, abs_grad)
