数字图像处理基础作业

合集下载

数字图像处理作业

1. Give a single intensity transformation function for spreading the intensities of an image so the lowest intensity is 0 and the highest is L-1. [为了扩展一幅图像的灰度，使其最低灰度为0、最高灰度为L-1，请给出一个单调的变换函数。

]Answer: Let f denote the original image. First subtract the minimum value of f denoted f min from f to yield a function whose minimum value is 0:Next divide g1 by its maximum value to yield a function in the range [0,1] and multiply the result by L一1 to yield a function with values in the range [0, L 一1]Keep in mind that f min is a scalar and f is an image.[让f表示原始图像。

首先从图像函数f中减掉f的最小值f min, 然后生成一个新的函数g1，它的最小值为0：接下来让g1的最大值除g1得到另一新的函数，它的值域在[0,1]区间，然后再乘上L一1，得到值域为[0, L一1]的新函数。

请注意f min是一个标量，而f是一个图像。

2.Explain why the discrete histogram equalization technique does not,in general,yield a flat histogram. [请解释为什么离散直方图均衡化技术一般不能得到平坦的直方图。

数字图像处理-作业题及部分答案解析

数字图像处理-作业题及部分答案解析1.数字图像与连续图像相比具有哪些优点?连续图像f(x,y与数字图像I(c,r中各量的含义是什么?它们有何联系和区别? (To be compared with an analog image, what are the advantages of a digital image? Let f(x,y be an analog image,I(r, c be a digital image, please give explanation and comparison for defined variables: f/I, x/r,and y/c2.图像处理可分为哪三个阶段? 它们是如何划分的?各有什么特点?(We can divide "imageprocessing"into 3 stages,what are they? how they are divided?What are their features?答:低级处理---低层操作,强调图像之间的变换,是一个从图像到图像的过程;中级处理---中层操作,主要对图像中感兴趣的目标进行检测和测量,从而建立对图像的描述,是一个从图像到数值或符号的过程;高级处理---高层操作,研究图像中各目标的性质和相互联系,得出对图像内容含义的理解以及对原来客观场景的解释;3.试从结构和功能等角度分析人类视觉中最基本的几个要素是什么?什么是马赫带效应? 什么是同时对比度?它们反映了什么共同问题? (According to the structure and function of the eyes,what are the basic elements in human vision? What is the Mach Band Effect? What is Simultaneous Contrast? What common facts can we infer from both Mach Band Effect and Simultaneous Contrast?答:人的视觉系统趋向于过高或过低估计不同亮度区域边界的现象称为“马赫带”效应;同时对比度指的是人的视觉系统对某个区域感觉到的亮度除了依赖于它本身的强度,还与背景有关.马赫带效应和同时对比度现象表明人所感觉到的亮度并不是强度的简单函数.4.比较说明像素邻域、连接、通路以及连通基本概念的联系与区别。

完整版数字图像处理作业题及部分答案

1.数字图像与连续图像相比具有哪些优点?连续图像f(x,y)与数字图像I(c,r)中各量的含义是什么?它们有何联系和区别? (To be compared with an analog image, what are the advantagesof a digital image? Let f(x,y) be an analog image, I(r, c) be a digital image, please giveexplanation and comparison for defined variables: f/I, x/r, and y/c)2.图像处理可分为哪三个阶段? 它们是如何划分的？各有什么特点? (We can divide image processing into 3 stages, what are they? how they are divided? What are their features?)答：低级处理---低层操作,强调图像之间的变换,是一个从图像到图像的过程;中级处理---中层操作,主要对图像中感兴趣的目标进行检测和测量,从而建立对图像的描述,是一个从图像到数值或符号的过程;高级处理---高层操作,研究图像中各目标的性质和相互联系,得出对图像内容含义的理解以及对原来客观场景的解释;3.试从结构和功能等角度分析人类视觉中最基本的几个要素是什么?什么是马赫带效应? 什么是同时对比度?它们反映了什么共同问题? (According to the structure and function of theeyes, what are the basic elements in human vision? What is the Mach Band Effect? What is Simultaneous Contrast? What common facts can we infer from both Mach Band Effect and Simultaneous Contrast?)答：人的视觉系统趋向于过高或过低估计不同亮度区域边界的现象称为“马赫带”效应;同时对比度指的是人的视觉系统对某个区域感觉到的亮度除了依赖于它本身的强度,还与背景有关. 马赫带效应和同时对比度现象表明人所感觉到的亮度并不是强度的简单函数.4.比较说明像素邻域、连接、通路以及连通基本概念的联系与区别。

硕士研究生《数字图像处理》作业

研究生《数字图像处理》考试1. 编写程序完成不同滤波器的图像频域降噪和边缘增强的算法并进行比较，得出结论。

● 图像频域降噪的实验原理与算法分析：图像的能量大部分集中在幅度谱的低频和中频部分，而图像的边缘和噪声对应于高频部分，因此能降低高频成分幅度的滤波器就能减弱噪声的影响，由卷积定理，在频域实现低通滤波的数学表达式：),(),(),(v u F v u H v u G =1. 理想低通滤波器（ILPF ）0),(),(01),(D v u D D v u D v u H >≤⎩⎨⎧=2. 巴特沃斯低通滤波器（BLPF ） nD v u D v u H 20),()12(11),(⎥⎦⎤⎢⎣⎡-+=3. 指数型低通滤波器（ELPF ） 2),(0),(nD v u D ev u H ⎥⎦⎤⎢⎣⎡-=● 图像频域降噪的实验过程： 1. 理想低通滤波器程序I=imread('xpy.jpg'); f=double(I); g=fft2(f); g=fftshift(g); [M,N]=size(g); d0=100;m=fix(M/2);n=fix(N/2); for i=1:Mfor j=1:Nd=sqrt((i-m)^2+(j-n)^2); if(d<=d0) h=1; else h=0; endresult(i,j)=h*g(i,j);endend>> result=ifftshift(result);>> J1=ifft2(result);>> J2=uint8(real(J1));>> imshow(J2)2.巴特沃斯低通滤波器程序I=imread('xpy.jpg');f=double(I);g=fft2(f);g=fftshift(g);[M,N]=size(g);nn=2;d0=30;m=fix(M/2);n=fix(N/2);for i=1:Mfor j=1:Nd=sqrt((i-m)^2+(j-n)^2);h=1/(1+0.414*(d/d0)^(2*nn));result(i,j)=h*g(i,j);endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));imshow(J2)3.高斯低通滤波器程序I=imread('xpy.jpg');f=double(I);g=fft2(f);g=fftshift(g);[M,N]=size(g);d0=100;m=fix(M/2);n=fix(N/2);for i=1:Mfor j=1:Nd=sqrt((i-m)^2+(j-n)^2);h=exp(-(d.^2)./(2*(d0^2)));result(i,j)=h*g(i,j);endendresult=ifftshift(result);J1=ifft2(result);J2=uint8(real(J1));imshow(J2)图像频域降噪的实验结果分析与讨论下面是理想低通滤波器、巴特沃斯低通滤波器、高斯低通滤波器的滤波效果分析与讨论。

《数字图像处理》作业题目

数字图像处理作业班级：Y100501姓名：**学号：*********一、编写程序完成不同滤波器的图像频域降噪和边缘增强的算法并进行比较，得出结论。

频域降噪。

对图像而言，噪声一般分布在高频区域，而图像真是信息主要集中在低频区，所以，图像降噪一般是利用低通滤波的方法来降噪。

边缘增强。

图像的边缘信息属于细节信息，主要由图像的高频部分决定，所以，边缘增强一般采取高通滤波，分离出高频部分后，再和原频谱进行融合操作，达到边缘增强，改善视觉效果，或者为进一步处理奠定基础的目的。

1频域降噪，主程序如下：I=imread('lena.bmp'); %读入原图像文件J=imnoise(I,'gaussian',0,0.02);%加入高斯白噪声A=ilpf(J,0.4);%理想低通滤波figure,subplot(222);imshow(J);title('加噪声后的图像');subplot(222);imshow(A);title('理想低通滤波');B=blpf(J,0.4,4);%巴特沃斯低通滤波subplot(223);imshow(B);title('巴特沃斯低通滤波');C=glpf(J,0.4);%高斯低通滤波subplot(224);imshow(C);title('高斯低通滤波');用到的滤波器函数的程序代码如下：function O=ilpf(J,p) %理想低通滤波，p是截止频率[f1,f2]=freqspace(size(J),'meshgrid');hd=ones(size(J));r=sqrt(f1.^2+f2.^2);hd(r>p)=0;y=fft2(double(J));y=fftshift(y);ya=y.*hd;ya=ifftshift(ya);ia=ifft2(ya);O=uint8(real(ia));function O=blpf(J,d,n) %巴特沃斯低通滤波器，d是截止频率，n是阶数[f1,f2]=freqspace(size(J),'meshgrid');hd=ones(size(J));r=f1.^2+f2.^2;for i=1:size(J,1)for j=1:size(J,2)t=r(i,j)/(d*d);hd(i,j)=1/(t^n+1);endendy=fft2(double(J));y=fftshift(y);ya=y.*hd;ya=ifftshift(ya);ia=ifft2(ya);O=uint8(real(ia));function O=glpf(J,D) %高斯滤波器，D是截止频率[f1,f2]=freqspace(size(J),'meshgrid');r=f1.^2+f2.^2;Hd=ones(size(J));for i=1:size(J,1)for j=1:size(J,2)t=r(i,j)/(D*D);Hd(i,j)=exp(-t);endendY=fft2(double(J));Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);ia=ifft2(Ya);O=uint8(real(ia));运行结果如图1所示。

数字图像处理作业题

数字图像处理作业题1. 以下是一个32级灰度图像（0表示黑色），其中包含了在灰色开始背景上的，带有一个白色标记的，直径是12英寸的黑色留声机唱片。

下面给出了图像的直方图。

试问象素间的间距是多大标记的尺寸是多大[0 0 0 0 100 200 2000 6000 2000 200 100 0 0 200 3000 9000 3000 200 0 0 50 100 400 100 50 0 0 0 0 0 0 0]解：像素数乘以像素间距的平方等于物体的面积（S ）。

表示唱片的像素总数：100+200+2000+6000+2000+200+100+200+50+100+400+100 +50=11300 S=22611300⨯=⨯πd d=(英寸) 表示白色标记的像素数为：50+100+400+100+50=700 S=22700r d ⨯=⨯π r=(英寸)2. 下面给出了在黑色背景上的白色台球的20级灰度图像的直方图0表示黑色），这个球是用每立方厘米克的材料制成的。

象素间距是1mm 。

试问球的重量是多少[0 100 500 3000 9000 3000 500 200 100 200 300 500 627 500 300 200 100 0 0 0]解：由直方图可知，表示台球的像素数为：100+200+300+500+627+500+300+200+100=2827S=222827r d ⨯=⨯π r=3cm 球的质量为： )(6.1695.1343g r M =⨯=π 原题：下面给出了在黑色背景上的白色台球的20级灰度图像的直方图0表示黑色），这个球是用每立方厘米克的材料制成的。

象素间距是1mm 。

试问球的重量是多少[0 200 500 3000 9000 3000 500 200 100 200 300 500 627 500 300 200 100 0 0 0]3.图像平滑的直观依据是什么不同的平滑方法是如何看待噪音并使用了何种改进以尽量降低其对边缘的模糊对于平滑的快速性和边缘保持，你有何见解解：图像在生成和传输过程中常受到各种噪声源的干扰和影响而使图像处理效果变差，反映在图像上，噪声使原本均匀和连续变化的灰度突然变大或减小，形成一些虚假的物体边缘或轮廓。

数字图像处理上机作业二

数字图像第二讲作业1.设计一个程序，对一幅灰度图像实现直方图均衡化处理。

画出均衡化前后的图像及其直方图.分析：要实现对一幅灰度图直方图的均衡化处理，须经四步来完成。

首先，获得直方图，及各灰度值对应像素数占总像素数的比例，用矩阵c表示，这是第一次直方图作业的内容。

第二步对c求累加和，得到s矩阵。

第三步，对s进行合理的近似，得到均衡图中的灰度等级。

这里我的处理办法是对s执行s*255操作，得到的值再取uint8(),这样就实现了对s*255的合理近似，即四舍五入，从而得到了新的灰度等级d。

第四步，赋予各像素点新的灰度值，计算出不同灰度等级对应的比例值Ps(sk)，这里我用：for j=1:256p(d(j)+1)=p(d(j)+1)+c(j);I2(find(I==j))=d(j);end来求得Ps(sk)，循环的作用是将灰度等级d(j)相同的各j对应的原灰度值比例c 累加起来。

I2(find(I==j))=d(j)是将灰度值为j的像素点的灰度值换为d(j);p(d(j)+1)中加1的作用是避免出现p(0)而进行的处理。

代码及注释如下：function junheng(x) % x为要分析的图像名加单引号I=imread(x);b=size(I);a=zeros(1,256); %a为一个1*256的矩阵分别记录灰度为0到255的像素的个数for m=1:b(1) %两个for语句将整张图的所有像素都扫描一遍for n=1:b(2)a(I(m,n)+1)= a(I(m,n)+1)+1; %将灰度为I(m,n)的像素个数存储在%a(I(m,n)+1)中,因为matlab里没有a(0)endendn=0:255;c=a/sum(a); %c为个灰度像素数占总像素数的比例subplot(2,2,1);bar(n,c);title('直方图Pr(rk)'); %画出直方图s=zeros(1,256);s(1)=c(1);for k=2:256s(k)=s(k-1)+c(k);end %s为Pr(rk)的累加和subplot(2,2,2);bar(n,s);title('sk'); %画出sk的图形d=uint8(s*255); %d是对s的合理近似值，即四舍五入p=zeros(1,256);I2=I*0;for j=1:256p(d(j)+1)=p(d(j)+1)+c(j);I2(find(I==j))=d(j);end %以d为灰度值，计算相同d值对应的c值的累加和f imwrite(I2,'junheng.bmp');subplot(2,2,3);bar(n,p);title('均衡化直方图Ps(sk)') %画出均衡化直方图figure;bar(n,c);title('直方图Pr(rk)'); %单独显示直方图figure;bar(n,p);title('均衡化直方图Ps(sk)') %单独显示均衡后直方图figure;subplot(1,2,1);imshow(x);title('原图'); %显示原图subplot(1,2,2);imshow('junheng.bmp');title('均衡后图');%显示均衡化之后图形运行：在命令窗口中输入junheng(‘Lenna.bmp’),输出如下四幅图：这幅图体现了整个程序设计的思路，先由直方图Pr(rk)到累加和的sk图形，最后到均衡图Ps(sk).这两幅图是单独显示的原直方图和均衡化后的直方图。

数字图像处理上机作业一.

数字图像处理上机作业一1.设计一个程序，绘制出一幅灰度图象的直方图。

Solution：代码及代码的说明：%作用：返回灰度矩阵a，并画出直方图function a=zhifangtu(x) % x为要分析的图像名加单引号I=imread(x);b=size(I);a=zeros(1,256); % a为一个1*256的矩阵分别记录灰度为0到255的像%的个数for m=1:b(1) %两个for语句将整张图的所有像素都扫描一遍for n=1:b(2)a(I(m,n)+1)= a(I(m,n)+1)+1; %将灰度为I(m,n)的像素个数存储在%a(I(m,n)+1)中,因为matlab里没有%a(0)endendn=0:255;bar(n,a);%画出直方图s=sum(a) %查看直方图的总的面积等于这张图的总像素值实验结果及分析：在命令窗口中输入zhifangtu('Lenna.bmp')返回s =262144，以及灰度矩阵a，同时有如下直方图输出：分析及结论：在命令窗口中用size命令可查知Lenna.bmp是512*512的，返回的s =262144恰等于512*512，说明所编的直方图的程序恰将所有的像素点都统计了，直方图的总面积等于像素总数。

直方图的作用也就是将一张图中不同灰度值对应像素数的一个统计。

在这个程序的编写中应注意a(I(m,n)+1)= a(I(m,n)+1)+1 不能写成a(I(m,n))= a(I(m,n))+1 ，应为在matlab中矩阵表示没有a(0),若某个像素点的灰度值是0，就会出错，故应写成a(I(m,n)+1)= a(I(m,n)+1)+1形式。

2.对同一场景但模糊程度不一样的三张数字图像绘制出其直方图, 计算每一幅图象所有像素灰度的方差。

图象的清晰度同灰度方差什么关系?Solution：代码及代码的说明：%作用：绘出模糊程度不一样的三张数字图像的直方图，并输出各自灰度方差I1=imread('tu1.bmp');I1=rgb2gray(I1); %转换为灰度图像imwrite(I1,'tu0.bmp'); %由于直方图只能对灰度图作用，故先将其转为灰度图subplot(2,2,1);zhifangtu('tu0.bmp');title('tu1直方图'); %绘出tu1.bmp的直方图k1=size(I1);I1=single(I1);I1=(I1-mean(mean(I1)')).^2; %个像素灰度值减去平均灰%值后再平方t1=sum(sum(I1)')/k1(1)/k1(2), %输出tu1.bmp的所有像素灰度的方差subplot(2,2,2);zhifangtu('tu2.bmp');title('tu2直方图'); %绘出tu2.bmp的直方图I2=imread('tu2.bmp');k2=size(I2);I2=single(I2);I2=(I2-mean(mean(I2)')).^2; %个像素灰度值减去平均灰度值后再%平方t2=sum(sum(I2)')/k2(1)/k2(2), %输出tu2.bmp的所有像素灰度的方差subplot(2,2,3);zhifangtu('tu3.bmp');title('tu3直方图'); %绘出tu3.bmp的直方图I3=imread('tu3.bmp');k3=size(I3);I3=single(I3);I3=(I3-mean(mean(I3)')).^2; %个像素灰度值减去平均灰度值%后再平方t3=sum(sum(I3)')/k3(1)/k3(2), %输出tu3.bmp的所有像素灰度的方差figure;subplot(2,2,1);imshow('tu1.bmp');title('tu1图'); %绘出tu1.bmp的图subplot(2,2,2);imshow('tu2.bmp');title('tu2图'); %绘出tu1.bmp的图subplot(2,2,3);imshow('tu3.bmp');title('tu3图'); %绘出tu1.bmp的图实验结果及分析：上述代码执行后，输出t1 = 7.3027e+003，t2= 6.5808e+003，t3=5.4860e+003；同时输出如下直方图：原始图：分析及结论：tu1，tu2，tu3三幅图是依次变模糊的，三张图的所有像素灰度方差依次为t1 = t1 = 7.3027e+003，t2= 6.5808e+003，t3=5.4860e+003，它们是依次变小的，可知图象的清晰度随灰度方差的变小而变得模糊。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

1.图像均值化2.图像的Gauss低通和Gauss高通3.对图像进行Gamma变化4.DCT变换，加上量化在反量化，和IDCT这四个题目，如果是对一些进行图像处理的程序员来讲或者很简单，但是我基本不接触图像处理这个方面的(虽然我头上挂着这个牌子)，基本原理不同，很难写程序。

不过幸好我同学是搞这个方面的，而且他的讲解能让我很快的知道我应该怎么去处理这个图像，而且通过网络的搜索，我发现wiki上的讲解真的相当的精准阿...不带让人纠结的数学公式，也不会长篇大论，有的就是通俗易懂的步骤和例子。

让人很快能知道我应该怎么都作就能完成这个效果的处理。

这是我一直很喜欢使用wiki的原因，推荐推荐阿....对于这个作业我本打算以最大的速度做完的，也不想真的去对原理进行真正的了解！所以只要知道怎么去做就可以了。

突然想到了那几天前我好像也学习QT,所以想连着这个一起作一次练习。

qt一个gui做的不错的库！而且简单的很...既然使用了QT那就要求我使用C++来写这个程序，C++这个语言，很久很久没使用了，主要是觉得自己对C++好像很陌生了，或者可以说是对面向对象这个思想的陌生。

关于这点，我也很想提出我的一点点想法。

虽然很早就开始使用面向对象去编程，可是好像自己一直没有入门面向对象这种思想。

使用C++的过程好像是将C++当成C来使用，很少很少说一定要使用到类，继承，多态这种特性。

看了很多书说，要学好C++就要放弃一些东西，把面向对象的这些东西学好。

可是面向对象到底是一种什么样的思想呢，有的时候很想用面向对象的思想去写程序，可是有时候却发现自己好像是为了面向对象而面向对象...好似纠结....应该是我经历的还不够吧。

关于这个作业，我最想讲的两个方面是：1.qt中的QImage这个对象，为了能让内存高效的访问，qt通过空间去换取时间的方法来提升效率。

让每一行都能被4B整除，这就是让qt本身会对每一行进行填充的过程，所以将一个一维的图像数据的转换为QImage是一个要小心的过程。

当然其实如果你懂的，那就是一个构造函数的问题，我在这方面犯过错，所以提出来了。

2.关于DCT的变换这个变换的公式很长很复杂，如果按公式来做的话，这个计算量是很大的，所以在网络搜索到一篇论文，好像是通过将其变为矩阵的乘法的运算，能减少相当大的运算量。

据说JPEG 也是用这种方式进行dct的变化的。

使用矩阵的运算简单很多，而且很多东西是很巧妙的联系着的。

关于dct应该是我花时间最多的一个方面，因为开始不懂是怎么样的过程，随便看了以下就开始写程序，最后不但程序有问题，连我改程序的心情都没有了。

直接把那个程序del了，重新写。

在写之前找了不少的论文和资料，在有很大把握了基础下，开始了编码的过程，那天花了5个小时才完成，最后还有一点点的暇眦，不过那个时候时间太晚了，坚持不住就去睡觉了。

第二天花了点时间调试以后，很顺利的完成了，而且变化的速度很快，比我使用公式法快很多。

而且效果很好的样子..让我很开心，这种成功的感觉很久没有了。

我只想把最好的那部分代码贴出来给大家看，DCT变换：当然如果大家需要这个运行的代码可以发email给我(smy19890720@),哈哈很多的bug我可不管的哦....[cpp]view plaincopyprint?1. void pimage::rgb_convert_yuv()2. {3. for(size_t j = 0; j < heigth8; j++)4. {5. for(size_t i = 0; i < width8; i++)6. {7. if(j < heigth && i < width)8. {9. int r = image[j*width*bpp+i*bpp + 2];10. int g = image[j*width*bpp+i*bpp + 1];11. int b = image[j*width*bpp+i*bpp];12. y[j][i] = 0.3 * r + g * 0.59 + b * 0.11;13. u[j][i] = -0.147 *r - 0.289 * g + 0.436 * b;14. v[j][i] = 0.615 * r - 0.515 * g - 0.100 * b;15. y[j][i] = y[j][i] - 128;16.17. //printf("r = %u,g = %u,b = %u\n",r,g,b);18. //printf("y = %d,u = %d,v = %d\n",y[j][i],u[j][i],v[j][i]);19. }20. else21. {22. y[j][i] = 0;23. u[j][i] = 0;24. v[j][i] = 0;25. }26. }27. }28. }29.30. v oid pimage::yuv_convert_rgb()31. {32. for(size_t i = 0; i < heigth8; i++)33. { for(size_t j = 0; j < width8; j++) {34. if(i < heigth && j < width)35. {36. int tmpy = y[i][j] + 128;37. int tmpu = u[i][j];38. int tmpv = v[i][j];39.40. int r = tmpy + (1.13983 * (tmpv));41. int g = tmpy - (0.58 * (tmpv)) - (0.39 *(tmpu));42. int b = tmpy + (2.03211 * (tmpu));43. if(r > 255)44. r = 255;45. if(g > 255)46. g = 255;47. if(b > 255)48. b = 255;49.50. if(r < 0)51. r = 0;52. if(g < 0)53. g = 0;54. if(b < 0)55. b = 0;56.57. //printf("r = %d,g = %d,b = %d\n",r,g,b);58. //printf("y = %d,u = %d,v = %d\n",tmpy,tmpu,tmpv);59.60.61. image[i*width*bpp+j*bpp] = b;62. image[i*width*bpp+j*bpp + 1] = g;63. image[i*width*bpp+j*bpp + 2] = r;64. }65. }66. }67. }68.69. v oid pimage::DCT()70. {71. int tmpy[8][8];72. int tmpu[8][8];73. int tmpv[8][8];74.75. rgb_convert_yuv();76. get_dct(dct,dct_t);77.78. for(int j = 0; j < heigth8/8; j++)79. {80. for(int i = 0; i < width8/8; i++)81. {82. for(int h = 0; h < 8; h++)83. {84. for(int w = 0; w < 8; w++)85. {86. tmpy[h][w] = y[j*8+h][i*8+w];87. tmpu[h][w] = u[j*8+h][i*8+w];88. tmpv[h][w] = v[j*8+h][i*8+w];89. //printf("j = %d,i = %d,h = %d,w = %d\n",j,i,h,w);90. }91. }92.93. matrix_mul(dct,tmpy,dct_t);94. matrix_mul(dct,tmpu,dct_t);95. matrix_mul(dct,tmpv,dct_t);96.97. lianhuay(tmpy);98. lianhuau(tmpu);99. lianhuau(tmpv);100. for(int h = 0; h < 8; h++) 101. {102. for(int w = 0; w < 8; w++) 103. {104. y[j*8+h][i*8+w] = tmpy[h][w]; 105. u[j*8+h][i*8+w] = tmpu[h][w]; 106. v[j*8+h][i*8+w] = tmpv[h][w]; 107. }108.109. }110. }111. }112. }113.114. void pimage::IDCT()115. {116. int tmpy[8][8];117. int tmpu[8][8];118. int tmpv[8][8];119.120.121. for(int j = 0; j < heigth8/8; j++) 122. {123. for(int i = 0; i < width8/8; i++) 124. {125. for(int h = 0; h < 8; h++) 126. {127. for(int w = 0; w < 8; w++) 128. {129. tmpy[h][w] = y[j*8+h][i*8+w]; 130. tmpu[h][w] = u[j*8+h][i*8+w]; 131. tmpv[h][w] = v[j*8+h][i*8+w]; 132. }133. }134.135. ilianhuay(tmpy);136. ilianhuau(tmpu);137. ilianhuau(tmpv);138.139. matrix_mul(dct_t,tmpy,dct);140. matrix_mul(dct_t,tmpu,dct);141. matrix_mul(dct_t,tmpv,dct);142.143. for(int h = 0; h < 8; h++)144. {145. for(int w = 0; w < 8; w++)146. {147. y[j*8+h][i*8+w] = tmpy[h][w];148. u[j*8+h][i*8+w] = tmpu[h][w];149. v[j*8+h][i*8+w] = tmpv[h][w];150. }151. }152.153. }154. }155. yuv_convert_rgb();156. }157. void pimage::get_dct(float (*dct)[8],float (*dct_T)[8]) 158. {159. for(int i = 0; i < 8; i++)160. {161. dct[0][i] = 1.0 / sqrt(8);162. dct_T[i][0] = dct[0][i];163. }164.165. for(int i = 1; i < 8; i++)166. for(int j = 0; j < 8; j++)167. {168. dct[i][j] = sqrt(2.0/8)*cos((2*j+1)*i*PI/(2.0*8));169. dct_T[j][i] = dct[i][j];170. }171. }172.173. void pimage::matrix_mul(float (*first)[8],int (*sec)[8],float (*third)[8]) 174. {175. float tmp[8][8];176. float tmp1[8][8];177.178. for(int i = 0;i < 8;i++)179. {180. for(int j = 0;j < 8;j++)181. {182. tmp[i][j] = 0;183. for(int k = 0;k < 8;k++)184. {185. tmp[i][j] += (first[i][k] * sec[k][j]);186. }187. }188. }189.190. for(int i = 0;i < 8; i++)191. {192. for(int j = 0;j < 8; j++)193. {194. tmp1[i][j] = 0;195. for(int k = 0; k < 8;k++)196. {197. tmp1[i][j] += (tmp[i][k] * third[k][j]);198. }199. }200. }201.202. for(int i = 0;i < 8;i++)203. {204. for(int j = 0;j < 8;j++)205. {206. sec[i][j] = round(tmp1[i][j]); 207. }208. }209. }210.211. int pimage::round(float thiz) 212. {213. if(fabs(thiz - (int)thiz) >= (1.0/2)) 214. {215. if(thiz > 0)216. {217. return (int)thiz + 1;218. }219. else220. {221. return (int)thiz - 1;222. }223. }224. else225. {226. return (int)thiz;227. }228. }229.230.231. void pimage::lianhuay(int (*thiz)[8]) 232. {233. static int matrix[8][8] = {234. {16,11,10,16,24,40,51,61}, 235. {12,12,14,19,26,58,60,55}, 236. {14,13,16,24,40,57,69,56}, 237. {14,17,22,29,51,87,80,62},238. {18,22,37,56,68,109,103,77},239. {24,35,55,64,81,104,113,92},240. {49,64,78,87,103,121,120,101},241. {72,92,95,98,112,100,103,99},242. };243.244. for(int i = 0; i < 8; i++)245. {246. for(int j = 0; j < 8; j++)247. {248. thiz[i][j] = round(thiz[i][j]*1.0/matrix[i][j]); 249. }250. }251. }252.253. void pimage::ilianhuay(int (*thiz)[8])254. {255. static int matrix[8][8] = {256. {16,11,10,16,24,40,51,61},257. {12,12,14,19,26,58,60,55},258. {14,13,16,24,40,57,69,56},259. {14,17,22,29,51,87,80,62},260. {18,22,37,56,68,109,103,77},261. {24,35,55,64,81,104,113,92},262. {49,64,78,87,103,121,120,101},263. {72,92,95,98,112,100,103,99},264. };265.266. for(int i = 0; i < 8; i++)267. {268. for(int j = 0; j < 8; j++)269. {270. thiz[i][j] = thiz[i][j] * matrix[i][j];271. }272. }273. }274.275. void pimage::lianhuau(int (*thiz)[8])276. {277. static int matrix[8][8] = {278. {17,18,24,47,99,99,99,99},279. {18,21,26,66,99,99,99,99},280. {24,26,56,99,99,99,99,99},281. {47,66,99,99,99,99,99,99},282. {99,99,99,99,99,99,99,99},283. {99,99,99,99,99,99,99,99},284. {99,99,99,99,99,99,99,99},285. {99,99,99,99,99,99,99,99},286. };287.288. for(int i = 0; i < 8; i++)289. {290. for(int j = 0; j < 8; j++)291. {292. thiz[i][j] = round(thiz[i][j]*1.0/matrix[i][j]); 293. }294. }295. }296.297. void pimage::ilianhuau(int (*thiz)[8])298. {299. static int matrix[8][8] = {300. {17,18,24,47,99,99,99,99},301. {18,21,26,66,99,99,99,99},302. {24,26,56,99,99,99,99,99},303. {47,66,99,99,99,99,99,99},304. {99,99,99,99,99,99,99,99},305. {99,99,99,99,99,99,99,99},306. {99,99,99,99,99,99,99,99},307. {99,99,99,99,99,99,99,99},308. };309.310. for(int i = 0; i < 8; i++)311. {312. for(int j = 0; j < 8; j++)313. {314. thiz[i][j] = thiz[i][j] * matrix[i][j];315. }316. }317. }318.319. void pimage::get_psnr()320. {321. uchar *tmp_image = (uchar*)calloc(width*heigth*bpp,sizeof(uchar)); 322. memcpy(tmp_image,image,width*heigth*bpp);323. DCT();324. IDCT();325.326. double psnr_r = 0;327. double psnr_g = 0;328. double psnr_b = 0;329.330. double mse_r = 0;331. double mse_g = 0;332. double mse_b = 0;333.334. for(int i = 0; i < this->heigth; i++)335. for(int j = 0; j < this->width; j++)336. {337. mse_b = mse_b + (image[i*this->width*bpp+j*bpp] - tmp_image[i* width*bpp+j*bpp]) * (image[i*this->width*bpp+j*bpp] - tmp_image[i*width*b pp+j*bpp]);338. mse_g = mse_g + (image[i*this->width*bpp+j*bpp + 1] - tmp_imag e[i*width*bpp+j*bpp + 1]) * (image[i*this->width*bpp+j*bpp + 1] - tmp_imag e[i*width*bpp+j*bpp + 1]);339. mse_r = mse_r + (image[i*this->width*bpp+j*bpp + 2] - tmp_image [i*width*bpp+j*bpp + 2]) * (image[i*this->width*bpp+j*bpp + 2] - tmp_image[i*width*bpp+j*bpp + 2]);340. }341. mse_r = mse_r / (heigth * width);342. mse_g = mse_g / (heigth * width);343. mse_b = mse_b / (heigth * width);344.345. psnr_r = 20 * log10(255/sqrt(mse_r));346. psnr_g = 20 * log10(255/sqrt(mse_g));347. psnr_b = 20 * log10(255/sqrt(mse_b));348.349. printf("r = %lf,g = %lf,b = %lf\n",psnr_r,psnr_g,psnr_b);350. free(tmp_image);351. }。