两种PSNR计算方法相差很大的一种可能原因——数据类型与计算精度

问题描述

有两个PSNR的计算代码。

代码1

def PSNR(img1, img2):
	pdb.set_trace()
	mse = np.mean( (img1 - img2) ** 2 )
	if mse == 0:
		return 100
	PIXEL_MAX = 255.0
	psnr = 20 * math.log10(PIXEL_MAX / math.sqrt(mse))
	print(psnr)
	return psnr

代码2

def getpsnr(img1,img2):
	pdb.set_trace()
	img1 = np.float64(img1) / (2**8-1)
	img2 = np.float64(img2) / (2**8-1)
	mse = np.mean(np.square(img1-img2))
	psnr = - 10 * np.log10(mse)
	print(psnr)
	return psnr

已知img1和img2的取值范围都是[0,255]. 显然，上面两种方法在数学上是等价的。然而，第一种方法计算出来的PSNR约为29dB，第二种方法的PSBR则为14dB左右，相差极大。.

 

问题分析

通过pdb分析具体的两张图片img1和img2，img1和img2的值如下：

(Pdb) img1
array([[[226, 195, 181],
        [231, 206, 193],
        [235, 207, 199],
        ...,
        [210, 179, 149],
        [211, 175, 148],
        [205, 171, 147]],

       [[234, 207, 193],
        [236, 209, 198],
        [237, 208, 199],
        ...,
        [209, 179, 151],
        [214, 179, 150],
        [212, 181, 147]],

       [[235, 209, 198],
        [236, 208, 199],
        [236, 209, 201],
        ...,
        [210, 181, 154],
        [213, 181, 151],
        [215, 182, 153]],

       ...,

       [[123, 131, 133],
        [119, 128, 130],
        [119, 128, 132],
        ...,
        [137, 142, 128],
        [136, 140, 126],
        [132, 135, 125]],

       [[125, 133, 132],
        [122, 129, 133],
        [120, 132, 134],
        ...,
        [139, 143, 127],
        [137, 140, 125],
        [132, 134, 122]],

       [[126, 135, 133],
        [123, 132, 135],
        [123, 132, 137],
        ...,
        [141, 144, 129],
        [138, 138, 129],
        [131, 133, 125]]], dtype=uint8)

(Pdb) img1.max()
254
(Pdb) img1.min()
13

(Pdb) img2
array([[[252, 239, 222],
        [245, 230, 211],
        [251, 230, 211],
        ...,
        [211, 174, 150],
        [246, 224, 200],
        [249, 228, 201]],

       [[250, 234, 216],
        [246, 229, 211],
        [250, 230, 211],
        ...,
        [213, 176, 151],
        [246, 226, 201],
        [247, 227, 199]],

       [[247, 229, 212],
        [246, 228, 208],
        [248, 230, 210],
        ...,
        [218, 181, 156],
        [248, 227, 203],
        [247, 227, 201]],

       ...,

       [[115, 114, 123],
        [115, 114, 123],
        [115, 114, 122],
        ...,
        [ 59,  66,  65],
        [ 67,  71,  73],
        [ 71,  75,  76]],

       [[116, 114, 125],
        [116, 114, 125],
        [115, 114, 123],
        ...,
        [ 69,  74,  74],
        [ 67,  71,  72],
        [ 72,  76,  77]],

       [[118, 116, 127],
        [117, 115, 126],
        [115, 113, 123],
        ...,
        [ 65,  69,  70],
        [ 68,  72,  73],
        [ 72,  76,  77]]], dtype=uint8)
(Pdb) img2.min()
0
(Pdb) img2.max()
255

第一种方法中的mse输出为：

(Pdb) mse = np.mean( (img1 - img2) ** 2 )
(Pdb) mse
69.70586649576823

第二种方法中的mse的输出为：

(Pdb) img11 = np.float64(img1) / (2**8-1)
(Pdb) img22 = np.float64(img2) / (2**8-1)
(Pdb) mse2 = np.mean(np.square(img11-img22))
(Pdb) mse2
0.013419700465140611

理论上，mse和mse2应该有如下关系：

mse2 = mse / (255*255)

但实际上：

(Pdb) mse / (255*255)
0.0010719856439180043
(Pdb) mse2
0.013419700465140611

另一方面，如果mse2 = mse / (255*255)，则方法2计算出的psnr = - 10 * np.log10(mse)将为：

# 若假设mse2 = mse / (255*255)，则方法2计算出的psnr将是

(Pdb) - 10 * np.log10(mse/(255*255))
29.69811030696365

# 方法1计算出的psnr为

(Pdb) 20 * math.log10(255 / math.sqrt(mse))
29.69811030696365

 

综上所述，两种算法的psnr结果有很大差异的原因是方法2中的mse2不满足

mse2 = mse / (255*255)

 

=============================================================================================

于是得到一个新的等价问题：

已知img1和img2是两个ndarray，为什么有：
(Pdb) mse1 = np.mean( (img1 - img2) ** 2 )
69.70586649576823
(Pdb) mse2 = np.mean( (img1/255 - img2/255) ** 2 )
0.013419700465140611
(Pdb) mse1/(255*255)
0.0010719856439180043 ≠ mse2

更进一步，发现

(Pdb) (img1**2)/((img1/1)**2)
array([[[2.58438406e-03, 3.60289283e-03, 7.60050060e-03],
        [2.11765147e-03, 4.61871995e-03, 3.46318022e-03],
        [3.34993210e-03, 2.26376345e-03, 4.46958410e-03],
        ...,
        [1.54195011e-03, 1.27961050e-03, 8.33295797e-03],
        [5.23348532e-03, 5.25714286e-03, 6.57414171e-03],
        [9.75609756e-04, 1.94931774e-03, 4.85908649e-03]],

       [[4.16392724e-03, 2.26376345e-03, 3.46318022e-03],
        [2.58546395e-03, 3.68581305e-03, 9.18273646e-04],
        [1.86935854e-03, 0.00000000e+00, 4.46958410e-03],
        ...,
        [3.68581305e-03, 1.27961050e-03, 7.45581334e-04],
        [4.97860075e-03, 1.27961050e-03, 1.01333333e-02],
        [3.20398718e-03, 7.60050060e-03, 4.85908649e-03]],

... ...

明显错误！

怀疑是数值精度问题，img1和img2的数据类型均为unit8

(Pdb) img1.dtype
dtype('uint8')

将方法1的代码改为

def PSNR(img1, img2):
	pdb.set_trace()
	img1 = np.float64(img1)
	img2 = np.float64(img2)
	mse = np.mean( (img1 - img2) ** 2 )
	if mse == 0:
		return 100
	PIXEL_MAX = 255.0
	psnr = 20 * math.log10(PIXEL_MAX / math.sqrt(mse))
	print(psnr)
	return psnr

发现方法1和方法2的计算结果一摸一样。

 

 

结论

方法1中img1和img2均为uint8类型，精度无法满足计算要求，故产生错误结果。

将方法1中的img1和img2改为np.float64后即可得到正确结果。

 

数值分析和数值精度并不是遥不可及的东西，一不小心就栽进去了，以后一定要小心，涉及到计算的最好早些归一化、早些转为float/couble类型。