如何在numpy.ndarray中计算x、y方向角度以得到预期结果？

场景与数组定义 #

创建了如下numpy.ndarray类型的surface数组：

import numpy as np

surface = np.zeros((10, 10))
surface[5, 5] = 1

数组内容为：

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]

预期输出 #

想要计算x方向的角度，预期输出为：

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 45 0 -45 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]

尝试的代码及结果 #

尝试了以下代码：

grad_y, grad_x = np.gradient(surface)
degrees_x = np.degrees(np.arctan(grad_x))

但得到的结果为：

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 26.560505 0 -26.560505 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]

解决方案 #

原因分析 #

np.gradient默认使用中心差分计算梯度，对于位置(5,4)，它的计算方式是(surface[5,5] - surface[5,3])/2 = (1-0)/2 = 0.5，arctan(0.5)对应的角度约为26.56度，和预期的45度不符——45度对应的是梯度值为1的情况（arctan(1)=45°）。

实现代码 #

要得到预期结果，需要计算不除以2的中心差分（即直接用右侧第二个元素减左侧第二个元素），代码如下：

import numpy as np

surface = np.zeros((10, 10))
surface[5, 5] = 1

# 初始化x方向梯度数组
grad_x = np.zeros_like(surface)
# 对中间列计算梯度：右侧第二个元素 - 左侧第二个元素
grad_x[:, 1:-1] = surface[:, 2:] - surface[:, :-2]
# 计算对应的角度
degrees_x = np.degrees(np.arctan(grad_x))

# 输出结果
print(degrees_x)

运行后，degrees_x的第5行即为[0. 0. 0. 0. 45. 0. -45. 0. 0. 0.]，完全匹配预期输出。

原理说明 #

• 对于位置(5,4)，梯度值为surface[5,5] - surface[5,3] = 1-0=1，arctan(1)转换为角度就是45°；
• 对于位置(5,6)，梯度值为surface[5,7] - surface[5,5] =0-1=-1，arctan(-1)转换为角度就是-45°；
• 其余位置梯度值为0，对应角度也为0，符合预期。

内容的提问来源于stack exchange，提问作者HJA24