OpenCV 区域编码和阈值编码实现图像压缩（8*8DCT变换，保留50%的系数）

《数字图像处理》的一次实验作业是要求如下：

我在网上找到了关于这个题目的Matlab版本的代码，链接如下：

http://wenku.baidu.com/link?url=6WlcVIE8sAg2Lnj6R7PQlv7tL7sNvBb7uhN3YZ37hj7xb855G_eA7Yu0ek34yl7UVKnRrlTUUBAOoNeoS5sb-wCMYwsfGgAoZHK-iLaU6rK

我参考上面的链接，主要做了下面两方面的改动：

（1）    将参考代码用OpenCV重新实现一遍（主要的难点在于OpenCV貌似没有类似于Matlab的blkproc那么方便的块处理函数，我采用的办法是定义一个8*8的感兴趣区域在图片

上面滑动）（我的代码在VS2013+OpenCV3.1上运行）；

（2）    参考代码中，阈值编码的部分有点莫名其妙，我采用的办法是：DCT变换之后对每个8*8的系数块取中值（系数取绝对值之后），绝对值小于中值的系数清零，每一个8*8的

块有不同的中值，也就会有不同的阈值。参考代码在阈值编码部分的关键代码如下：

上面的程序用8*8的矩阵g小于中值的系数置为0之后作为掩膜，不太理解为什么要用g矩阵，即使是这样，c=median(b)这句程序也有问题，举个例子，比如：

b=[-3 -2 -1 0 1 2 3]，那么median(b)为0，显然达不到去除小于中值的系数的效果，所以应该把c=median(b)改成c=median(abs(b))

我的程序的主要思路如下：

（1）    以灰度图方式读入原始图片，并将图片转换成double类型；

（2）    用自行编写的blkproc_DCT函数，对原始图像进行8*8分块的DCT变换；

（3）    用自行编写的regionalCoding和thresholdCoding函数，分别对DCT变换后的系数矩阵进行区域编码和阈值编码；

（4）    用自行编写的IDCT函数，从regionalCoding和thresholdCoding函数得到的结果中恢复出原始图片，显示并比较结果。

下面是具体的代码（运行的时候只需要修改一下imread后面的图片路径即可）：

#include<iostream>
#include<opencv2/opencv.hpp>

using namespace std;
using namespace cv;

void blkproc_DCT(Mat);	//功能等同于Matlab的blkproc函数（blkproc(I,[8 8],'P1*x*P2',g,g')），用于对图像做8*8分块DCT
Mat blkproc_IDCT(Mat);	//功能等同于Matlab的blkproc函数（blkproc(I,[8 8],'P1*x*P2',g',g)），用于做8*8分块DCT逆变换，恢复原始图像
void regionalCoding(Mat);	//区域编码函数，功能等同于Matlab的blkproc函数（blkproc(I1,[8 8],'P1.*x',a)）
void thresholdCoding(Mat);	//阈值编码函数，功能等同于Matlab的blkproc函数（blkproc(I1,[8 8],'P1.*x',a)）
double get_medianNum(Mat &);	//获取矩阵的中值，用于阈值编码

#define M_PI 3.141592653

int main()
{
	Mat ucharImg = imread("F:\\My_Desktop\\data\\frame\\hand_test\\hand_test_0.jpg",0);	//以灰度图的形式读入原始的图像
	imshow("srcImg", ucharImg);
	Mat doubleImg;
	ucharImg.convertTo(doubleImg, CV_64F);	//将原始图像转换成double类型的图像，方便后面的8*8分块DCT变换
	blkproc_DCT(doubleImg);	//对原图片做8*8分块DCT变换
	//分别进行区域编码和阈值编码
	Mat doubleImgRegion, doubleImgThreshold;
	doubleImgRegion = doubleImg.clone();
	doubleImgThreshold = doubleImg.clone();
	regionalCoding(doubleImgRegion);	//对DCT变换后的系数进行区域编码
	thresholdCoding(doubleImgThreshold);	//对DCT变换后的系数进行阈值编码
	//进行逆DCT变换
	Mat ucharImgRegion, ucharImgThreshold;
	ucharImgRegion = blkproc_IDCT(doubleImgRegion);	
	imshow("RegionalCoding", ucharImgRegion);
	ucharImgThreshold = blkproc_IDCT(doubleImgThreshold);
	imshow("ThresholdCoding", ucharImgThreshold);
	waitKey(0);
}

void blkproc_DCT(Mat doubleImgTmp)
{
	Mat ucharImgTmp;
	Mat DCTMat = Mat(8, 8, CV_64FC1);	//用于DCT变换的8*8的矩阵
	Mat DCTMatT;	//DCTMat矩阵的转置
	Mat ROIMat = Mat(8, 8, CV_64FC1);	//用于分块处理的时候在原图像上面移动
	double a = 0, q;	//DCT变换的系数
	for (int i = 0; i < DCTMat.rows; i++)
	{
		for (int j = 0; j < DCTMat.cols; j++)
		{
			if (i == 0)
			{
				a = pow(1.0 / DCTMat.rows, 0.5);
			}
			else
			{
				a = pow(2.0 / DCTMat.rows, 0.5);
			}
			q = ((2 * j + 1)*i*M_PI) / (2 * DCTMat.rows);
			DCTMat.at<double>(i, j) = a*cos(q);
		}
	}
	DCTMatT = DCTMat.t();
	//ROIMat在doubleImgTmp以8为步长移动，达到与Matlab中的分块处理函数blkproc相同的效果
	//此程序中，若图片的高或者宽不是8的整数倍的话，最后的不足8的部分不进行处理
	int rNum = doubleImgTmp.rows / 8;
	int cNum = doubleImgTmp.cols / 8;
	for (int i = 0; i < rNum; i++)
	{
		for (int j = 0; j < cNum; j++)
		{
			ROIMat = doubleImgTmp(Rect(j * 8, i * 8, 8, 8));
			ROIMat = DCTMat*ROIMat*DCTMatT;
		}
	}
	doubleImgTmp.convertTo(ucharImgTmp, CV_8U);
	imshow("DCTImg", ucharImgTmp);
}

Mat blkproc_IDCT(Mat doubleImgTmp)
{
	//与blkproc_DCT几乎一样，唯一的差别在于：ROIMat = DCTMatT*ROIMat*DCTMat（转置矩阵DCTMatT和DCTMat交换了位置）
	Mat ucharImgTmp;
	Mat DCTMat = Mat(8, 8, CV_64FC1);	
	Mat DCTMatT;	
	Mat ROIMat = Mat(8, 8, CV_64FC1);	
	double a = 0, q;	
	for (int i = 0; i < DCTMat.rows; i++)
	{
		for (int j = 0; j < DCTMat.cols; j++)
		{
			if (i == 0)
			{
				a = pow(1.0 / DCTMat.rows, 0.5);
			}
			else
			{
				a = pow(2.0 / DCTMat.rows, 0.5);
			}
			q = ((2 * j + 1)*i*M_PI) / (2 * DCTMat.rows);
			DCTMat.at<double>(i, j) = a*cos(q);
		}
	}
	DCTMatT = DCTMat.t();
	int rNum = doubleImgTmp.rows / 8;
	int cNum = doubleImgTmp.cols / 8;
	for (int i = 0; i < rNum; i++)
	{
		for (int j = 0; j < cNum; j++)
		{
			ROIMat = doubleImgTmp(Rect(j * 8, i * 8, 8, 8));
			ROIMat = DCTMatT*ROIMat*DCTMat;
		}
	}
	doubleImgTmp.convertTo(ucharImgTmp, CV_8U);
	return ucharImgTmp;
}

void regionalCoding(Mat doubleImgTmp)
{
	int rNum = doubleImgTmp.rows / 8;
	int cNum = doubleImgTmp.cols / 8;
	Mat ucharImgTmp;
	Mat ROIMat = Mat(8, 8, CV_64FC1);	//用于分块处理的时候在原图像上面移动
	for (int i = 0; i < rNum; i++)
	{
		for (int j = 0; j < cNum; j++)
		{
			ROIMat = doubleImgTmp(Rect(j * 8, i * 8, 8, 8));
			for (int r = 0; r < ROIMat.rows; r++)
			{
				for (int c = 0; c < ROIMat.cols; c++)
				{
					//8*8块中，后四行置0
					if (r>4)
					{
						ROIMat.at<double>(r, c) = 0.0;
					}
				}
			}
		}
	}
	doubleImgTmp.convertTo(ucharImgTmp, CV_8U);
	imshow("regionalCodingImg", ucharImgTmp);
}

void thresholdCoding(Mat doubleImgTmp)
{
	int rNum = doubleImgTmp.rows / 8;
	int cNum = doubleImgTmp.cols / 8;
	double medianNumTmp = 0;
	Mat ucharImgTmp;
	Mat ROIMat = Mat(8, 8, CV_64FC1);	//用于分块处理的时候在原图像上面移动
	for (int i = 0; i < rNum; i++)
	{
		for (int j = 0; j < cNum; j++)
		{
			ROIMat = doubleImgTmp(Rect(j * 8, i * 8, 8, 8));
			medianNumTmp = get_medianNum(ROIMat);
			for (int r = 0; r < ROIMat.rows; r++)
			{
				for (int c = 0; c < ROIMat.cols; c++)
				{
					if (abs(ROIMat.at<double>(r,c))<0)
					{
						ROIMat.at<double>(r, c) = 0;
					}
				}
			}
		}
	}
	doubleImgTmp.convertTo(ucharImgTmp, CV_8U);
	imshow("thresholdCodingImg", ucharImgTmp);
}

double get_medianNum(Mat & imageROI)	//获取矩阵的中值
{
	vector<double> vectorTemp;
	double tmpPixelValue = 0;
	for (int i = 0; i < imageROI.rows; i++)	//将感兴趣区域矩阵拉成一个向量
	{
		for (int j = 0; j < imageROI.cols; j++)
		{
			vectorTemp.push_back(abs(imageROI.at<double>(i, j)));
		}
	}
	for (int i = 0; i < vectorTemp.size() / 2; i++)	//进行排序
	{
		for (int j = i + 1; j < vectorTemp.size(); j++)
		{
			if (vectorTemp.at(i) > vectorTemp.at(j))
			{
				double temp;
				temp = vectorTemp.at(i);
				vectorTemp.at(i) = vectorTemp.at(j);
				vectorTemp.at(j) = temp;
			}
		}
	}
	return vectorTemp.at(vectorTemp.size() / 2 - 1);	//返回中值
}