数字图像处理100问—37 PSNR

提示：内容整理自：https://github.com/gzr2017/ImageProcessing100Wen
CV小白从0开始学数字图像处理

37 PSNR

离散余弦逆变换中如果不使用8作为系数，而是使用4作为系数的话，图像的画质会变差。来求输入图像和经过离散余弦逆变换之后的图像的峰值信噪比吧！再求出离散余弦逆变换的比特率

峰值信噪比（Peak Signal to Noise Ratio）缩写为PSNR，用来表示信号最大可能功率和影响它的表示精度的破坏性噪声功率的比值，可以显示图像画质损失的程度。

峰值信噪比越大，表示画质损失越小。峰值信噪比通过下式定义。MAX表示图像点颜色的最大数值。如果取值范围是[0,255]的话，那么MAX的值就为255。MSE表示均方误差（Mean Squared Error），用来表示两个图像各个像素点之间差值平方和的平均数。

PSNR = 10 * log10(MAX^2 / MSE)
MSE = Sum_{y=0:H-1} Sum_{x=0:W-1} (I1(x,y) - I2(x,y))^2 / (HW)

如果我们进行8x8的离散余弦变换，离散余弦逆变换的系数为KxK的话，比特率按下式定义：

bitrate = 8 * K^2 / 8^2

代码如下：

1.引入库

CV2计算机视觉库

import cv2
import numpy as np
import matplotlib.pyplot as plt

2.读入数据

img = cv2.imread("imori.jpg").astype(np.float32)
H, W, C = img.shape

3.灰度化

gray = 0.2126 * img[..., 2] + 0.7152 * img[..., 1] + 0.0722 * img[..., 0]

4.DCT

T = 8
K = 4
X = np.zeros((H, W), dtype=np.float64)
#indx = np.tile(np.arange(T), (T, 1))
#indy = np.arange(T).repeat(T).reshape(T, -1)
#dct = np.ones_like(indx, dtype=np.float32)
#dct[:, 0] /= np.sqrt(2)
#dct[0] /= np.sqrt(2)

def w(x, y, u, v):
    cu = 1.
    cv = 1.
    if u == 0:
        cu /= np.sqrt(2)
    if v == 0:
        cv /= np.sqrt(2)
    theta = np.pi / (2 * T)
    return (( 2 * cu * cv / T) * np.cos((2*x+1)*u*theta) * np.cos((2*y+1)*v*theta))
    
for yi in range(0, H, T):
    for xi in range(0, W, T):
        for v in range(T):
            for u in range(T):
                for y in range(T):
                    for x in range(T):
                        X[v+yi, u+xi] += gray[y+yi, x+xi] * w(x,y,u,v)
                """
                _x = indx + xi * T
                _y = indy + yi * T
                _u = u + xi * T
                _v = v + yi * T
                X[_v, _u] = np.sum(C * gray[_y, _x] * np.cos((2*indx+1)*u*np.pi/(2*T)) * np.cos((2*indy+1)*v*np.pi/(2*T)))
                """

5.IDCT

out = np.zeros((H, W), dtype=np.float64)

for yi in range(0, H, T):
    for xi in range(0, W, T):
        for y in range(T):
            for x in range(T):
                for v in range(K):
                    for u in range(K):
                        out[y+yi, x+xi] += X[v+yi, u+xi] * w(x,y,u,v)
                """
                _u = indx + xi * T
                _v = indy + yi * T
                _x = x + yi * T
                _y = y + xi * T
                out[_y, _x] = np.sum(C * X[_v, _u] * np.cos((2*x+1)*indx*np.pi/(2*T)) * np.cos((2*y+1)*indy*np.pi/(2*T))) * 4. / (T ** 2)
                """

out[out>255] = 255
out = np.floor(out).astype(np.uint8)

# MSE
v_max = 255.
mse = np.sum(np.power(np.abs(gray.astype(np.float32) - out.astype(np.float32)), 2)) / (H * W)
psnr = 10 * np.log10(v_max ** 2 / mse)

print("PSNR >>", psnr)

bitrate = 1. * T * K ** 2 / (T ** 2)
print("bitrate >>", bitrate)

6.保存结果

cv2.imshow("result", out)
cv2.waitKey(0)
cv2.imwrite("out.jpg", out)

7.结果