虽然Python是弱类型编程语言,但是不同的数据类型对计算结果有很大的影响。这点和Matlab不同。
下面以numpy.gradient()函数对一幅图像做梯度运算来说明。
python源码如下,是Feature extraction & image processing for computer vision (3rd edition) 中 code 4.17 的实现:
import cv2
import numpy as np
import matplotlib.pyplot as plt
def GradCorner(img,op='tan'):
# img = np.mat(img,dtype=np.float32)
row,col = img.shape
img_out = np.zeros((row,col))
Mx, My= np.gradient(img,edge_order=2)
# Mx = cv2.Sobel(img,cv2.CV_32FC1,1,0)
# My = cv2.Sobel(img,cv2.CV_32FC1,0,1)
EdgeX = cv2.Sobel(img,cv2.CV_32FC1,1,0)
EdgeY = cv2.Sobel(img,cv2.CV_32FC1,0,1)
Mag = np.sqrt(EdgeX**2+EdgeY**2)
# plt.subplot(144);plt.imshow(Mag);plt.title('Edge');plt.show()
A = np.arctan2(-Mx,My)
# max supression ,to be impleted...
Mxx,Mxy = np.gradient(Mx,1)
Myx,Myy = np.gradient(My,1)
# comput curvature
for x in np.arange(col):
for y in np.arange(row):
if Mag[y,x]!=0:
My2=My[y,x]**2;Mx2=Mx[y,x]**2;MxMy=Mx[y,x]*My[y,x]
if Mx2+My2!=0:
if op=='Ti':
img_out[y,x]=(1/(Mx2+My2)**1.5)*(My2*Mxx[y,x]-MxMy*Myx[y,x]
-Mx2*Myy[y,x]+MxMy*Mxy[y,x])
elif op=='N':
img_out[y, x] = (1 / (Mx2 + My2) ** 1.5) * (Mx2 * Myx[y, x]-MxMy * Mxx[y, x]
- MxMy * Myy[y, x] + My2 * Mxy[y, x])
elif op=='NI':
img_out[y, x] = (1 / (Mx2 + My2) ** 1.5) *(-Mx2 * Myx[y, x]+MxMy * Mxx[y, x]
- MxMy * Myy[y, x] + My2 * Mxy[y, x])
else :
img_out[y, x] = (1 / (Mx2 + My2) ** 1.5) * (-My2 * Mxx[y, x] - MxMy * Myx[y, x]
+ Mx2 * Myy[y, x] - MxMy * Mxy[y, x])
plt.subplot(221);plt.imshow(img);plt.title('origin')
plt.subplot(222);plt.imshow(Mx);plt.title('grad X')
plt.subplot(223);plt.imshow(My);plt.title('grad Y')
plt.subplot(224);plt.imshow(img_out);plt.title('Corner');plt.show()
# plt.tight_layout()
return img_out
if __name__ == '__main__':
image = cv2.imread('u.tif',flags=0)
image = cv2.resize(image,(100,100))
corimg = GradCorner(image,op='TI')
# cv2.imshow('orign',image)
# cv2.imshow('gradcorner',corimg)
# cv2.waitKey()
plt.imshow(corimg)
plt.show()
如果第10行是采用numpy 的gradient函数:
Mx, My= np.gradient(img,edge_order=2)
那么求梯度的结果是:
显然求出的梯度和最终的角点都不正确:不论是X方向还是Y方向,从255(红)变0(蓝)的梯度都没有计算出来。查找gradient实现的源码,有这样一段代码(function_base.py, 1077-1084行):
# Use second order differences where possible
out = np.empty_like(y, dtype=otype)
slice1[axis] = slice(1, -1)
slice2[axis] = slice(2, None)
slice3[axis] = slice(None, -2)
# 1D equivalent -- out[1:-1] = (y[2:] - y[:-2])/2.0
out[slice1] = (y[slice2] - y[slice3])/2.0
其实,out的计算就是将原来的输入矩阵y 在某一个维度偏移了两个位置,求其差,作为梯度。这一步后,out的计算就不正确了。Debug时发现out是uint8类型,差值中的负值-255会被算为1。查看out的元素,也印证了这一想法。如果把最后一句改为:
out[slice1] = (y[slice3] - y[slice2])/2.0
那么没有出现的梯度就出现了,而原来出现的又没有了。
所以在自己的代码中加了:
img = np.mat(img,dtype=np.float32)
将uint8类型强制转换为float32,最终的结果就正确了:
在图像处理中,cv2.imread的返回结果大多是unsigned类型,虽然在python代码中没有直接指出,但是在运算的时候,python还是严格遵循了不同类型运算的法则。
Python是弱类型,但不是没有数据类型。
这一点需要引起大家的注意。