基于块变换和哈尔小波编码的

1、医学图像处理实验报告实验十一：基于块变换和哈尔小波编码的图像压缩与重建摘要 本次实验的实验目的及主要内容是：Ø 基于块变换的图像压缩重建Ø 基于小波变换的图像压缩重建一、技术讨论1.1实验原理图像信息熵：它表征图像灰度分布的聚集特性，却不能反映图像灰度分布的空间特征，为了表征这种空间特征，可以在一维熵的基础上引入能够反映灰度分布空间特征的特征量来组成图像的二维熵。一阶熵：相邻两个像素的像素值和作为统计标准；二阶熵：相邻三个像素的像素值和作为统计标准。保真度准则：表征电子设备输出再现输入信号的相似程度。保真度越高，无线电接收机输出的声音或电视机输出的影像越逼真。块变换编码:映

2、射变换：函数变换，常用的称为正交变换-量化-编码小波编码：小波变换（不损失数据，它是EZW算法具有渐进性的基础）-量化（对数据产生损失，损失程度取决于量化阈值的大小）-熵编码（对每个输入字精确的确定它的概率，并根据概率生成一个和合适的代码，使输出的码流小雨输入码流）。1.2实验方法1.void perform( const Mat &InputArray, Mat &OutputArray,const int reserve,const int flag);method->perform(originalimage,dst,1,DCT); /depend on dct o

3、r dft2. 调用对象的成员函数 Imagentropy来求信息熵，ZEROORDER为零阶熵，FIRST_ORDER为一阶熵，method->Imagentropy(originalimage,FIRST_ORDER);3. 使用对象method的成员函数FidelityCriteria求保真度；二、结果与讨论2.1实验结果1. 改变块变换编码中所使用变换的种类（DFT与DCT），结果图如下：DFT： 图一（a） （b）DCT： 图二（a） （b）2. 改变哈尔小波编码中分级参数level与置零系数阈值res，结果图如下：Level=1 res=1 图三（a） （b） Level=1

4、 res=0.01 图四（a） （b）Level=2 res=0.01 图五（a） （b）Level=3 res=0.01 图六（a） （b）3.了解：可利用块变换编码的程序计算图像的0、1、2阶信息熵，通过改变主程序中计算信息熵的函数最后一个输入参数ZERO_ORDER为FIRST_ORDER或SECOND_ORDER实现。2.2实验讨论1. 块变换编码中，比较使用不同变换运行程序后出现的均方根误差ems与信噪比SNR的结果值，简单解释出现差别的原因，并对两种变换的压缩效果做出评价；由图一和图二可以看出，DCT变换结果优于DFT。离散傅里叶变换DFT是傅里叶变换在时域和频域上都呈离散的形式。

5、实际应用中通常采用快速傅里叶变换计算DFT。离散余弦变换DCT是与傅里叶变换相关的一种变换，DCT相当于一个长度大概是它两倍的DFT，但DCT是对一个实偶函数进行变换的。从形式上来看，DCT是线性可逆变换。当图像压缩比增大时，压缩效率减小，图像质量降低。但是，当压缩比较高时，这种块的独立处理带来了压缩效率上的限制和方块效应问题，方块效应成为限制DCT编码质量的主要因素。2.试阐释哈尔小波编码中分级参数level与置零系数阈值res对结果的影响；根据图四、图五、图六的对照试验可以看出：分级参数对于图像变换后的清晰度有很大影响，level值越大，则图像越模糊。这和理论分析结果相互对应；根据图三和图

6、四的对照试验可以看出：res的值对图像的压缩效果不明显（由于实验代码问题），理论分析出置零系数阈值是对图像压缩影响最大的系数。小于res的值置零，大于res的值不变，因此res值越大，压缩效果越明显。附录（实验代码）.pro文件程序：#-# Project created by QtCreator 2014-04-04T14:35:39#-QT += coreQT -= guiTARGET = 2dfourieCONFIG += consoleCONFIG -= app_bundleTEMPLATE = appSOURCES += main.cppINCLUDEPATH+=d:Qtopencv

7、2.2includeopencvd:Qtopencv2.2includeopencv2d:Qtopencv2.2includeLIBS+=d:Qtopencv2.2liblibopencv_calib3d220.dll.ad:Qtopencv2.2liblibopencv_contrib220.dll.ad:Qtopencv2.2liblibopencv_core220.dll.ad:Qtopencv2.2liblibopencv_features2d220.dll.ad:Qtopencv2.2liblibopencv_flann220.dll.ad:Qtopencv2.2liblibopen

8、cv_gpu220.dll.ad:Qtopencv2.2liblibopencv_highgui220.dll.ad:Qtopencv2.2liblibopencv_imgproc220.dll.ad:Qtopencv2.2liblibopencv_legacy220.dll.ad:Qtopencv2.2liblibopencv_ml220.dll.ad:Qtopencv2.2liblibopencv_objdetect220.dll.ad:Qtopencv2.2liblibopencv_video220.dll.aWavelet文件程序如下：头文件code.h:#ifndef CODE_H#

9、define CODE_H#include <iostream>#include <opencv2/core/core.hpp>#include <opencv2/highgui/highgui.hpp>#include <opencv2/imgproc/imgproc.hpp>using namespace cv;Mat WDT( const Mat &_src, const String _wname, const int _level );void wavelet( const String _wname, Mat &_lo

10、wFilter, Mat &_highFilter );Mat IWDT( const Mat &_src, const String _wname, const int _level );Mat waveletDecompose( const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter );Mat waveletReconstruct( const Mat &_src, const Mat &_lowFilter, const Mat &_highFilter

11、 );#endif/ CODE_HCode.cpp:#include"code.h"/ 小波变换Mat WDT( const Mat &_src, const String _wname, const int _level ) /int reValue = THID_ERR_NONE; Mat src = Mat_<float>(_src); Mat dst = Mat:zeros( src.rows, src.cols, src.type() ); int N = src.rows; int D = src.cols; / 高通低通滤波器 Mat lo

12、wFilter; Mat highFilter; wavelet( _wname, lowFilter, highFilter ); / 小波变换 int t=1; int row = N; int col = D; while( t<=_level ) /先进行行小波变换 for( int i=0; i<row; i+ ) / 取出src中要处理的数据的一行 Mat oneRow = Mat:zeros( 1,col, src.type() ); for ( int j=0; j<col; j+ ) oneRow.at<float>(0,j) = src.at&

13、lt;float>(i,j); oneRow = waveletDecompose( oneRow, lowFilter,highFilter ); / 将src这一行置为oneRow中的数据 for ( int j=0; j<col; j+ ) dst.at<float>(i,j) = oneRow.at<float>(0,j); / 小波列变换 for ( intj=0; j<col; j+) / 取出src数据的一列输入 Mat oneCol = Mat:zeros( row, 1,src.type() ); for ( int i=0; i&l

14、t;row; i+ ) oneCol.at<float>(i,0)= dst.at<float>(i,j); oneCol=( waveletDecompose( oneCol.t(),lowFilter,highFilter ).t(); for ( inti=0; i<row; i+ ) dst.at<float>(i,j) = oneCol.at<float>(i,0); / 更新 row /= 2; col /=2; t+; src=dst; normalize(dst,dst,0,1,CV_MINMAX); return dst;

15、/ 小波逆变换Mat IWDT( const Mat &_src, const String _wname, const int _level ) /int reValue = THID_ERR_NONE; Mat src= Mat_<float>(_src); Mat dst=Mat:zeros(src.rows,src.cols,src.type() ); int N =src.rows; int D = src.cols; /高通低通滤波器 Mat lowFilter; Mat highFilter; wavelet( _wname,lowFilter, highFi

16、lter ); /小波变换 intt=1; int row =N/std:pow( 2., _level-1); int col= D/std:pow(2., _level-1); while ( row<=N&&col<=D ) /小波列逆变换 for ( int j=0;j<col;j+) /取出src数据的一行输入 Mat oneCol =Mat:zeros( row,1,src.type() ); for ( int i=0; i<row;i+) oneCol.at<float>(i,0)=src.at<float>(i,

17、j); oneCol = ( waveletReconstruct( oneCol.t(),lowFilter,highFilter ) ).t(); for(int i=0;i<row;i+) dst.at<float>(i,j) = oneCol.at<float>(i,0); /行小波逆变换 for( int i=0;i<row; i+ ) /取出src中要处理的数据的一行 Mat oneRow = Mat:zeros( 1,col,src.type(); for ( intj=0; j<col; j+ ) oneRow.at<float&

18、gt;(0,j) =dst.at<float>(i,j); oneRow=waveletReconstruct(oneRow, lowFilter, highFilter ); /将src这一行置为oneRow中的数据 for( int j=0;j<col; j+ ) dst.at<float>(i,j) =oneRow.at<float>(0,j); row*=2; col *= 2; src=dst; return dst;/*.*/ 调用函数/ 生成不同类型的小波，现在只有haar，sym2void wavelet(constString _wn

19、ame, Mat &_lowFilter, Mat &_highFilter) if( _wname="haar"|_wname="db1" ) intN=2; _lowFilter = Mat:zeros(1, N,CV_32F ); _highFilter =Mat:zeros(1, N, CV_32F ); _lowFilter.at<float>(0,0)=1/sqrtf(N); _lowFilter.at<float>(0,1)=1/sqrtf(N); _highFilter.at<float>

20、;(0,0)=-1/sqrtf(N); _highFilter.at<float>(0, 1)= 1/sqrtf(N); if (_wname ="sym2" ) int N=4; floath =-0.483,0.836, -0.224, -0.129; floatl=-0.129,0.224, 0.837,0.483; _lowFilter =Mat:zeros( 1,N,CV_32F ); _highFilter=Mat:zeros(1, N,CV_32F ); for (int i=0;i<N;i+) _lowFilter.at<float&

21、gt;(0,i)= li; _highFilter.at<float>(0,i)=hi; /小波分解Mat waveletDecompose(const Mat &_src, const Mat &_lowFilter, constMat &_highFilter ) assert(_src.rows=1 &&_lowFilter.rows=1&&_highFilter.rows=1 ); assert(_src.cols>=_lowFilter.cols &&_src.cols>=_highFil

22、ter.cols ); Mat src = Mat_<float>(_src); intD= src.cols; Mat lowFilter=Mat_<float>(_lowFilter); Mat highFilter=Mat_<float>(_highFilter); / 频域滤波，或时域卷积；ifft(fft(x) *fft(filter)=cov(x,filter) Mat dst1=Mat:zeros( 1, D, src.type(); Mat dst2 =Mat:zeros( 1,D, src.type() ); filter2D(src,ds

23、t1,-1,lowFilter ); filter2D(src,dst2,-1, highFilter); / 下采样 Mat downDst1 =Mat:zeros(1,D/2,src.type(); Mat downDst2 =Mat:zeros(1,D/2, src.type() ); resize(dst1, downDst1,downDst1.size(); resize(dst2, downDst2, downDst2.size() ); / 数据拼接 for ( inti=0;i<D/2; i+) src.at<float>(0, i)=downDst1.at&

24、lt;float>(0, i); src.at<float>(0, i+D/2)=downDst2.at<float>(0, i); return src;/小波重建Mat waveletReconstruct(const Mat &_src,constMat &_lowFilter, constMat &_highFilter) assert( _src.rows=1 &&_lowFilter.rows=1&& _highFilter.rows=1 ); assert(_src.cols>=_lowF

25、ilter.cols && _src.cols>=_highFilter.cols ); Mat src=Mat_<float>(_src); int D =src.cols; Mat lowFilter = Mat_<float>(_lowFilter); Mat highFilter =Mat_<float>(_highFilter); /插值; Mat Up1 = Mat:zeros(1,D,src.type() ); Mat Up2=Mat:zeros(1, D,src.type();/ /插值为0/ for( inti=0,c

26、nt=1;i<D/2; i+,cnt+=2)/ / Up1.at<float>(0, cnt)=src.at<float>(0,i); /<前一半/ Up2.at<float>(0, cnt) =src.at<float>(0, i+D/2 );/< 后一半/ /线性插值 Mat roi1( src, Rect(0,0,D/2,1); Mat roi2( src, Rect(D/2,0,D/2, 1); resize(roi1,Up1, Up1.size(),0,0, INTER_CUBIC ); resize(roi2,Up2

27、,Up2.size(), 0,0,INTER_CUBIC ); / 前一半低通，后一半高通 Mat dst1=Mat:zeros(1, D,src.type() ); Mat dst2= Mat:zeros(1,D, src.type() ); filter2D(Up1,dst1,-1, lowFilter); filter2D( Up2, dst2, -1, highFilter ); /结果相加 dst1 =dst1+ dst2; normalize(dst1,dst1,0,1,CV_MINMAX); returndst1;Main.cpp:#include"code.h&quo

28、t;intmain(int argc, char *argv) int level=1; Mat src=imread("d:/lena.tif",CV_LOAD_IMAGE_GRAYSCALE); Mat dst=WDT(src,"haar",level); float res=0.01; Mat te = Mat_<float>(dst); for (int i=round(dst.rows/(pow(2,level);i<dst.rows;i+) for (int j=round(dst.cols/(pow(2,level);j&

29、lt;dst.cols;j+) te.at<float>(i,j)=(abs(te.at<float>(i,j)<res)*te.at<float>(i,j); Mat th=IWDT(te,"haar",level); imshow("dst",dst); imshow("th",th); waitKey(0); return 0;Transform-coding文件程序如下：头文件transformcoding.h#ifndefTRANFORM_CODING_H#define TRANFO

30、RM_CODING_H#include <iostream>#include <opencv2/core/core.hpp>#include <opencv2/highgui/highgui.hpp>#include <opencv2/imgproc/imgproc.hpp>using namespace cv;using namespace std;enum DCT=0,DFT=1,ZERO_ORDER=2,FIRST_ORDER=3,SECOND_ORDER=4;class Tranformcodingpublic: double ems,S

31、NR,H; /构造函数，初始化成员变量 Tranformcoding() this->ems=0; this->SNR=0; this->H=0; this->index=0; Tranformcoding() delete this; void FidelityCriteria(const Mat &compressed, const Mat &decompressed); void perform( const Mat &InputArray, Mat &OutputArray,const int reserve,const int

32、flag); void Imagentropy(const Mat &InputArray,const int flag);private: void Cosinecoding(Mat &InputArray,Mat &Z,const int reserve); void Fouriercoding( Mat &InputArray, Mat &Z,const int reserve); int index; Point maxidex;#endif/TRANFORM_CODING_H源文件main.cpp:#include "tranform

33、coding.h"int main(int argc, char *argv) Mat originalimage = imread("/home/timlentse/Qtproject/lena.tif",CV_LOAD_IMAGE_GRAYSCALE); / Tranformcoding是我定义的一个类，使用这个类来实例化一个对象,名字叫method Tranformcoding *method=new Tranformcoding; Mat dst; /前面我们已经有了method这个对象了，当然可以使用它来调用它本身的成员函数，来执行我们预先定义好的行为

34、/例如method->Imagentropy 执行信息熵的计算 method->Imagentropy(originalimage,ZERO_ORDER); / i want to output the second_order entroy cout<<"-the entropy ="<<method->H<<endl; / method->perform 执行块转换编码 method->perform(originalimage,dst,8,DFT); /depend on dct or dft imsh

35、ow("the original image",originalimage); imshow ("the effect", dst); waitKey(0); /method->FidelityCriteria 执行保真度计算 method->FidelityCriteria(originalimage,dst); cout<<"-the ems ="<<method->ems<<endl<<"-and the SNR ="<<metho

36、d->SNR<<endl; return 0;Transformcoding.cpp：#include "tranformcoding.h"/这是计算图像保真度（包括ems,SNR)的子程序，你需要输入两幅图像，即原图和压缩的图/第一个参数是原图，第二个参数是压缩后的图void Tranformcoding:FidelityCriteria(const Mat &compressed, const Mat &decompressed) if (compressed.size != decompressed.size| compressed.

37、type()!= decompressed.type() cout<<"-"<<endl; cout<<"error:nothing done because the two images are not in the same size or type "<<endl <<"It is recommended to check your two images"<<endl; cout<<"-"<<endl; retur

38、n ; SNR = norm(decompressed, NORM_L2); ems = norm(compressed, decompressed, NORM_L2); SNR = pow(SNR/ems), 2); ems = sqrt(ems*ems/(compressed.rows* compressed.cols);/以下是求图像信息熵的子函数 ,flag标志位：/-ZERO ORDER ：求零阶熵/-FIRST ORDER ：求一阶熵/-SECOND ORDER ：求二阶熵void Tranformcoding:Imagentropy(const Mat &InputArr

39、ay, const int flag) int array766= 0; /像素值统计部分，放在array数组里面 int row= InputArray.rows; int col= InputArray.cols; if (InputArray.isContinuous() col= col*row; row= 1; for (int i(0); i<row; +i) for(intj(0); j<col; +j) const uchar*p =InputArray.ptr<uchar>(i); if (flag=2) index=pj; else if(flag=

40、 3) index=pj+(j!=col-1)*pj+1; else if(flag= 4) index=(j!=0)*pj-1+pj+(j!=col-1)*pj+1; arrayindex+; /计算信息熵 double pr(0); for(intk(0); k<766; k+) pr= arrayk; if (pr!=0) pr/=(InputArray.cols* InputArray.rows); H=- (log(pr)/log(2)*pr + H; /定义perform函数来响应用户的需求，perform函数是用来实施块转换编码的/根据标志位flag和保留系数个数来执行块转

41、换编码voidTranformcoding:perform( constMat &InputArray, Mat &OutputArray, const int reserve, const int flag) OutputArray= Mat_<float>(InputArray); float t88= 16,11,10,16,24,40,51,61, 12,12,14,19,26,58,60,55, 14,13,16,24,40,57,69,56, 14,17,22,29,51,87,80,62, 18,22,37,56,68,109,103,77, 24,3

42、5,55,64,81,104,113,92, 49,64,78,87,103,121,120,101, 72,92,95,98,112,100,103,99; Mat Z = Mat(8,8,CV_32F,t); for(intu(0); u<OutputArray.rows/8;+u) for (int v(0);v<OutputArray.cols/8; +v) Mat temp(OutputArray,Rect(v*8,u*8,8,8); if (flag=1) Fouriercoding(temp,Z, reserve); elseif(flag=0) Cosinecoding(temp, Z, reserve); else std:cout<<"-no matching fuctioncall,maybethe flagiswrong-"<<std:endl; OutputArray.conv