Python 数据类型对运算的影响

虽然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是弱类型，但不是没有数据类型。
这一点需要引起大家的注意。