求输入图像和经过离散余弦逆变换之后的图像的峰值信噪比。并求出离散余弦逆变换的比特率。
一、名词简介
DCT- 离散余弦变换，在（声音、图像）数据压缩中得到了广泛的使用。
PSNR-峰值信噪比（Peak Signal to Noise Ratio）缩写为PSNR，用来表示信号最大可能功率和影响它的表示精度的破坏性噪声功率的比值，可以显示图像画质损失的程度。峰值信噪比越大，表示画质损失越小。
MSE-均方误差，（mean-square error, MSE）是反映估计量与被估计量之间差异程度的一种度量。
二、参考c++代码
#include <opencv2/core.hpp>#include <opencv2/highgui.hpp>#include <iostream>#include <math.h>#include <complex>const int height = 128, width = 128, channel = 3;// DCT超参数 DCT hyper-parameterint T = 8;int K = 4;// DCT coefficientstruct dct_str {double coef[height][width][channel];};// 离散余弦变换 Discrete Cosine transformationdct_str dct(cv::Mat img, dct_str dct_s){double I;double F;double Cu, Cv;for(int ys = 0; ys < height; ys += T){for(int xs = 0; xs < width; xs += T){for(int c = 0; c < channel; c++){for(int v = 0; v < T; v ++){for(int u = 0; u < T; u ++){F = 0;if (u == 0){Cu = 1. / sqrt(2);} else{Cu = 1;}if (v == 0){Cv = 1. / sqrt(2);}else {Cv = 1;}for (int y = 0; y < T; y++){for(int x = 0; x < T; x++){I = (double)img.at<cv::Vec3b>(ys + y, xs + x)[c];F += 2. / T * Cu * Cv * I * cos((2. * x + 1) * u * M_PI / 2. / T) * cos((2. * y + 1) * v * M_PI / 2. / T);}}dct_s.coef[ys + v][xs + u][c] = F;}}}}}return dct_s;}// 逆离散余弦变换 Inverse Discrete Cosine transformationcv::Mat idct(cv::Mat out, dct_str dct_s){double f;double Cu, Cv;for (int ys = 0; ys < height; ys += T){for (int xs = 0; xs < width; xs += T){for(int c = 0; c < channel; c++){for (int y = 0; y < T; y++){for (int x = 0; x < T; x++){f = 0;for (int v = 0; v < K; v++){for (int u = 0; u < K; u++){if (u == 0){Cu = 1. / sqrt(2);} else {Cu = 1;}if (v == 0){Cv = 1. / sqrt(2);} else { Cv = 1;}f += 2. / T * Cu * Cv * dct_s.coef[ys + v][xs + u][c] * cos((2. * x + 1) * u * M_PI / 2. / T) * cos((2. * y + 1) * v * M_PI / 2. / T);}}f = fmin(fmax(f, 0), 255);out.at<cv::Vec3b>(ys + y, xs + x)[c] = (uchar)f;}}}}}return out;}// 计算均方误差 Compute MSEdouble MSE(cv::Mat img1, cv::Mat img2){double mse = 0;for(int y = 0; y < height; y++){for(int x = 0; x < width; x++){for(int c = 0; c < channel; c++){mse += pow(((double)img1.at<cv::Vec3b>(y, x)[c] - (double)img2.at<cv::Vec3b>(y, x)[c]), 2);}}}mse /= (height * width);return mse;}// 计算峰值信噪比 Compute PSNRdouble PSNR(double mse, double v_max){return 10 * log10(v_max * v_max / mse);}// 计算比例 Compute bitratedouble BITRATE(){return T * K * K / T * T;}// Mainint main(int argc, const char* argv[]){double mse;double psnr;double bitrate;// read original imagecv::Mat img = cv::imread("imori.jpg", cv::IMREAD_COLOR);// DCT coefficientdct_str dct_s;// output imagecv::Mat out = cv::Mat::zeros(height, width, CV_8UC3);// DCTdct_s = dct(img, dct_s);// IDCTout = idct(out, dct_s);// MSE, PSNRmse = MSE(img, out);psnr = PSNR(mse, 255);bitrate = BITRATE();std::cout << "MSE: " << mse << std::endl;std::cout << "PSNR: " << psnr << std::endl;std::cout << "bitrate: " << bitrate << std::endl;cv::imwrite("out.jpg", out);//cv::imshow("answer", out);//cv::waitKey(0);cv::destroyAllWindows();return 0;}
三、参考python代码
import cv2import numpy as npimport matplotlib.pyplot as plt# DCT hyoer-parameterT = 8K = 4channel = 3# DCT weightdef 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))# DCTdef dct(img):H, W, _ = img.shapeF = np.zeros((H, W, channel), dtype=np.float32)for c in range(channel):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):F[v+yi, u+xi, c] += img[y+yi, x+xi, c] * w(x,y,u,v)return F# IDCTdef idct(F):H, W, _ = F.shapeout = np.zeros((H, W, channel), dtype=np.float32)for c in range(channel):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, c] += F[v+yi, u+xi, c] * w(x,y,u,v)out = np.clip(out, 0, 255)out = np.round(out).astype(np.uint8)return out# MSEdef MSE(img1, img2):H, W, _ = img1.shapemse = np.sum((img1 - img2) ** 2) / (H * W * channel)return mse# PSNRdef PSNR(mse, vmax=255):return 10 * np.log10(vmax * vmax / mse)# bitratedef BITRATE():return 1. * T * K * K / T / T# Read imageimg = cv2.imread("imori.jpg").astype(np.float32)# DCTF = dct(img)# IDCTout = idct(F)# MSEmse = MSE(img, out)# PSNRpsnr = PSNR(mse)# bitratebitrate = BITRATE()print("MSE:", mse)print("PSNR:", psnr)print("bitrate:", bitrate)# Save resultcv2.imshow("result", out)cv2.waitKey(0)cv2.imwrite("out.jpg", out)
四、OpencvSharp版本部分代码
// 读取原始文件Mat img = Cv2.ImRead(@"C:\Users\xiaomao\Desktop\lena.jpg", ImreadModes.Grayscale);// 转为32FC1Mat src = new Mat();img.ConvertTo(src, MatType.CV_32FC1);// 进行离散余弦变换Mat dct = new Mat(img.Size(), MatType.CV_32FC1);Cv2.Dct(src, dct);// 进行逆离散余弦变换Mat idct = new Mat(img.Size(), MatType.CV_32FC1);Cv2.Idct(dct, idct);// 将逆离散余弦变换结果转换为8UC1idct.ConvertTo(src, MatType.CV_8UC1);// 计算psnr，使用lena女神的图像测试结果为361.9909921956double psnr = Cv2.PSNR(img, idct);// 存储变换后的图像Cv2.ImWrite(@"C:\Users\zyh\Desktop\lena_idct.jpg", src);
