离散余弦函数（c++实现）

理论：

图像处理中常用的正交变换除了傅里叶变换外，还有其他一些有用的正交变换,其中离散余弦就是一种。离散余弦变换表示为DCT( Discrete Cosine Transformation)，常用于图像处理和图像识别等。

一维离散余弦变换

正变换


(1)


(2)

式中F(u)是第u个余弦变换系数，u是广义频率变量,u=1,2,3……N-1; f(x)是时域N点序列, x=0,1,2……N-1

反变换


（3）

显然，式(1)式(2)和式(3)构成了一维离散余弦变换对。

二维离散余弦变换

正变换


（4）

式(4)是正变换公式。其中f(x,y)是空间域二维向量之元素， x,y=0,1,2,……N-1；F(u,v)是变换系数阵列之元素。式中表示的阵列为N×N

反变换


（5）

式中的符号意义同正变换式一样。式(4)和式(5)是离散余弦变换的解析式定义。

矩阵表示法

更为简洁的定义方法是采用矩阵式定义。根据以上公式定义可知，离散余弦变换的系数矩阵可以写成如下：



如果令N＝4，那么由一维解析式定义可得如下展开式。



写成矩阵式



若定义F(u)为变换矩阵，A为变换系数矩阵，f(x)为时域数据矩阵，则一维离散余弦变换的矩阵定义式可写成如下形式

[F(u)]=[A][f(x)] （6）

同理，可得到反变换展开式



写成矩阵式即

[f(x)]=[A]T[F(u)] （7）

二维离散余弦变换也可以写成矩阵式：

[F(u,v)]=[A][f(x,y)][A]T （8）

[f(x,y)]=[A]T[F(u,v)][A]

式中[f(x,y)]是空间数据阵列，A是变换系数阵列，[F(u,v)]是变换矩阵，[A]T是[A]的转置。

对二维图像进行离散余弦变换

由以上对二维离散余弦变换的定义及公式（7）可知，求二维图像的离散余弦变换要进行以下步骤：

1.获得图像的二维数据矩阵f(x,y)；

2.求离散余弦变换的系数矩阵[A]；

3.求系数矩阵对应的转置矩阵[A]T；

4.根据公式（7）[F(u,v)]=[A][f(x,y)][A]T 计算离散余弦变换；

实现代码：

#include<opencv2/opencv.hpp>
#include<bitset>
using namespace std;
using namespace cv;

void coefficient(const int& n,double** p,double** q){
double t = 1.0 / sqrt(n+0.0);
for (int i = 0; i < n;i++){
q[0][i] = t;
p[i][0] = t;
}
for (int i = 1; i < n;i++){
for (int j = 1; j < n;j++){
q[i][j] = sqrt(2.0 / n)*cos(i*(j + 0.5)*CV_PI / n);
p[j][i] = q[i][j];
}
}
}
void matrix_computation(double** result,double** q,double** p, int &n){
double t = 0;
for (int i = 0; i < n;i++){
for (int j = 0; j < n;j++){
t = 0;
for (int k = 0; k < n;k++){
t += q[i][j] * p[j][i];
}
result[i][j] = t;
}
}
}
void DCT(Mat_<uchar> img,double** iMatrix,int &n){
for (int i = 0; i < n;i++){
for (int j = 0; j < n; j++){
iMatrix[i][j] = img(i, j);
}
}

double** coefficient1 = new double*
;
double** coefficient2 = new double*
;
double** tmp = new double*
;
for (int i = 0; i < n; i++){
coefficient1[i] = new double
;
coefficient2[i] = new double
;
tmp[i] = new double
;
}
coefficient(n, coefficient1, coefficient2);
matrix_computation(tmp,coefficient1,iMatrix,n);
matrix_computation(iMatrix, coefficient2, tmp, n);

for (int i = 0; i < n; i++){
delete[] coefficient1[i];
delete[] coefficient2[i];
delete[] tmp[i];
}
delete[] coefficient1;
delete[] coefficient2;
delete[] tmp;
}
void calc(double** iMatrix,bitset<64> &it){
double sum=0;
for (int i = 0; i < 8; i++){
for (int j = 0; j < 8; j++){
sum += iMatrix[i][j];

}
}
double average = sum / 64;
for (int i = 0; i < 8; i++){
int pos = i * 8;
for (int j = 0; j < 8; j++){
it.at(pos+j)=iMatrix[i][j] >= average ? 1 : 0;
}
}
}
int main(){
Mat img = imread("6.jpg");
bitset<64> bits;

resize(img, img, Size(32, 32));
cvtColor(img, img, CV_BGR2GRAY);
int n = img.rows;
double** imatrix = new double*
;
for (int i = 0; i < n;i++){
imatrix[i] = new double
;
}
DCT(img,imatrix,n);
calc(imatrix, bits);
for (int i = 0; i < 64; i++){
cout << bits.at(i);
}
cout << endl;
Mat img1 = imread("1.jpg");
bitset<64> bits1;
resize(img1, img1, Size(32, 32));
cvtColor(img1, img1, CV_BGR2GRAY);
int n1 = img1.rows;
double** imatrix1 = new double*[n1];
for (int i = 0; i < n1; i++){
imatrix1[i] = new double[n1];
}
DCT(img1, imatrix1, n1);
calc(imatrix1, bits1);
for (int i = 0; i < 64; i++){
cout << bits1.at(i);
}

for (int i = 0; i < n1; i++){
delete[] imatrix1[i];
}

for (int i = 0; i < n; i++){
delete[] imatrix[i];
}
delete[] imatrix;
cin.get();
delete[] imatrix1;
return 0;
}