Python Numpy gradient源码解析

复习图像梯度，发现Numpy有一个梯度计算函数，解析它的源码和需要注意的问题，最后自定义一个梯度函数

---

目录

1. 用法解析
2. 示例和问题
3. 源码解析
4. 自定义

---

用法解析

参考：numpy.gradient

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修改计算方式

差分公式

参考：Finite difference

一阶前向差分：

$$△f(x_k)=f(x_{k+1})-f(x_k)$$

一阶后向差分：

$$▽f(x_k)=f(x_k)-f(x_{k-1})$$

一阶中心差分：

$$\delta f(x_k)=\frac{f(x_{k+1})-f(x_{k-1})}{2}$$

二阶中心差分：

$$f^{\prime \prime }(x_k)=f^{\prime }(x_k)-f^{\prime }(x_{k-1})=f(x_{k+1})-f(x_k)-(f(x_k)-f(x_{k-1}))=f(x_{k+1})+f(x_{k-1})-2\cdot f(x_k)$$

---

示例和问题

示例和问题

计算一维或二维数组的梯度，以及计算图像梯度

一维数组

>>> 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，得到负数会溢出，所以需要先转换成更高精度的类型

参考：Data types

>>> 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. ]])

---

源码解析

np.gradient函数源码：gradient source

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，计算图像梯度如下：

参考：python matplotlib 显示图像

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)