基于奇异值分解的图像压缩与除噪

一、本报告所用的一些基本原理

第一，虑噪过程。

巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。在振幅的对数对角频率的波得图上，从某一边界角频率开始，振幅随着角频率的增加而逐步减少，趋向负无穷大。

本程序通过不断调整通带截止频率、阻带下限截止频率、通带衰减及阻带衰减，力图使滤噪后的图像逼近未加噪图像。

 

第二，图像压缩思路。

图像信息的数据量非常的大，随着各种成像设备的分辨率不断提高，单幅图像所包含的数据量也越来越大，大数据量的图像信息会给存储器的存储容量、通信信道的带宽以及计算机的处理速度增加极大的压力。为了解决这个问题，必须对图像进行压缩处理

在本次数值分析的图像压缩课程模拟实践当中，原图是一个2000*3000的一个矩阵，这里就是一个10^6的一个很大的数据量，我们就要想办法在保证图像质量的前提下，尽可能的压缩图像，在实践当中分别选择了1/10、1/100、1/1000、1/10000进行练习。

我们小组在压缩实践中主要运用的思路是基于SVD算法的图像压缩，在进行1/10的压缩时，我们对2000*3000的矩阵进行20*30的分块，每一小块为一个100*100的矩阵，进行奇异值分解时，我们只选取了前5个奇异值，这样在传输过程中的数据量仅为K=100*5+5*5+5*100=1025，K/（100*100）近似等于1/10，这样在图像的传输过程中，就只需对U（:,1:5)、S(1:5,1:5)、V(:,1:5)'三个矩阵的数据进行传输，而不必在传输以前的矩阵数据，数据接收方只需对这三个矩阵进行处理，即可得到压缩了1/10后的图像。

在进行1/100的压缩时，对滤波后的图像矩阵整体进行奇异值分解，这一次取前12个奇异值，这样的数据量K=2000*12+12*12+12*3000=60144，K/(2000*3000)近似于1/100。

这样的思路在进行1/1000及更高比例的压缩时遇到了问题，因为要进行1/1000的压缩，只能取1个奇异值左后才能够满足压缩要求，但是这样压缩之后的图像却严重的失真了，原因当然是清楚的，对于2000个奇异值只取最大的1个奇异值，图像的很多东西就被丢弃了，这样图像压缩遇到了瓶颈，不得不停滞下来思考其他方法；经过小组成员之间的讨论和与其他小组的交流中，有了新的思路，先将2000*3000的滤波后矩阵进行200*300的分块，每一小块有100个数据，通过c=mean（CB（:))函数对每一小块的矩阵进行均值化处理，这样一个10*10的矩阵变成了一个均值c，再通过新矩阵赋值CD(i,j)=c,最后一个2000*3000的矩阵通过均值化处理再赋值变成了200*300的矩阵，在保持了图像的基本质量中图像已经压缩了1/100，后面只需对这200*300的矩阵压缩1/10即可达到压缩1/1000的要求，在显示回图像的过程中，因为2000*3000的矩阵变成了200*300的矩阵，因此x、y坐标的长度的缩小十分之一（这个问题一直没有注意到，导致压缩程序老是出现错误，显示不出图像，后再小组成员的讨论中终于寻找出了错误）

最后是进行1/10000的压缩，也就是对上一步中的那个200*300的矩阵压缩1/100，进行奇异值分解，只取前2个奇异值，这样K=200*2+2*2+2*300=1004，K/（200*300）近似于1/60，已和1/100处在同一个数量级，能够满足要求。

 

二、关于图像噪点过滤、压缩及误差所编写程序

clear all ;

clc

 

%% 生成2000*3000数据

f = @(x,y) exp(x-2*x.^2-y.^2).*sin(6*(x + y +x.*y.^2));

x =linspace(-1,1,3e3);

y =linspace(-1.2,1.2,2e3); y = y(:) ;

A = f(repmat(x,2e3,1),repmat(y,1,3e3)) ;

AA=A;

r = -0.2 + (0.2+0.2).*rand(2000,3000); %增加噪声项

A=A+r;

 

%% 原始数据图

figure,mesh(x,y,A); %三维网格图

xlabel('原始加噪点图');

view([-40,80]);     %调整视角

 

 

%% Butterworth对图像进行滤波

disp('１．Butterworth Filtering');

%低通滤波器技术要求，经上分析，拟定

%通带截止频率为1.5Hz，阻带下限截止频率为4Hz

%通带衰减为0.7dB，阻带衰减为1.3dB

fs=3000;Wp=2*pi*1.5/fs;Ws=2*pi*4/fs;Rp=0.7;Rs=1.3;

Omip=Wp/pi;Omis=Ws/pi;                 %归一化技术要求

[N,Wn]=buttord(Omip,Omis,Rp,Rs);       %确定滤波器的阶数

disp(['The order of Butterworth Filtering is',num2str(N)]);

[b,a]=butter(N,Wn);                 %确定Butterworth滤波器转移函数系数向量

R=filter(b,a,A);

figure,mesh(x,y,R); %三维网格图,滤波后图像

xlabel('滤波后的图像');

view([-40,80]);     %调整视角

VC=AA(:);  %将矩阵变成列矢量

VC0=R(:); %将矩阵变成列矢量

RMSE0=sqrt(sum((VC-VC0).^2)/6000000); %求均方根误差

 

 

%% 分块奇异值压缩1/10

BB=R;

for i=1:1:20

    for j=1:1:30

       BC=(BB(((i-1)*100+1):i*100,((j-1)*100+1):j*100));

        [U,S,V]=svd(BC);                 %对每个100*100块做奇异值分解

       BD=U(:,1:5)*S(1:5,1:5)*V(:,1:5)';  %取前5个特征值

       BK(((i-1)*100+1):i*100,((j-1)*100+1):j*100)=BD;

    end

end

 

%% 分块奇异值分解后的图

figure,mesh(x,y,BK);

xlabel('压缩1/10后的图像');

view([-40,80]);

VC1=BK(:); %将矩阵变成列矢量

RMSE1=sqrt(sum((VC-VC1).^2)/6000000); %求均方根误差

 

 

%% 奇异值分解进行压缩1/100

[U,S,V]=svd(R);     %奇异值分解出U,S,V三个矩阵

C = U(:,1:12)*S(1:12,1:12)*V(:,1:12)'; %取U的前12列，S的前12个奇异值，V的前12列

U2= U(:,1:12);

S2= S(1:12,1:12);

V2= V(:,1:12)';

 

%% 奇异值分解后的图

figure,mesh(x,y,C);

xlabel('压缩1/100后的图像');

view([-40,80]);

VC2=C(:); %将矩阵变成列矢量

RMSE2=sqrt(sum((VC-VC2).^2)/6000000); %求均方根误差

 

 

%% 分块奇异值10*10均值后压缩1/1000

CA=R;

for i=1:1:200

    for j=1:1:300

       CB=(CA(((i-1)*10+1):i*10,((j-1)*10+1):j*10));

        c=mean(CB(:));%全部平均

                CD(i,j)=c;

    end

end

 

%% 分块奇异值压缩1/10

CE=CD;

for i=1:1:2

    for j=1:1:3

       CF=(CE(((i-1)*100+1):i*100,((j-1)*100+1):j*100));

        [U,S,V]=svd(CF);                 %对每个10*10块做奇异值分解

       CG=U(:,1:5)*S(1:5,1:5)*V(:,1:5)';  %取前5个特征值

       CH(((i-1)*100+1):i*100,((j-1)*100+1):j*100)=CG;

    end

end

 

%% 分块奇异值压缩后的图

x =linspace(-1,1,3e2);

y =linspace(-1.2,1.2,2e2); y = y(:) ;

figure,mesh(x,y,CH);

xlabel('压缩1/1000后的图像');

view([-40,80]);

 

%% 压缩后的CD奇异值分解进行压缩1/10000

E=CD;

[U,S,V]=svd(E);    %奇异值分解出U,S,V三个矩阵

D = U(:,1:2)*S(1:2,1:2)*V(:,1:2)';  %取U的前2列，S的前2个奇异值，V的前2列

U4= U(:,1:2);

S4= S(1:2,1:2);

V4= V(:,1:2)';

 

%% 奇异值分解后的图

figure,mesh(x,y,D);

xlabel('压缩1/10000后的图像');

view([-40,80]);

 

 

三、程序运行后的结果见附件二

评价：对图像进行压缩，不可避免的要引入失真。我们要做的就是在图像信号的最终用户察觉不出或能够忍受这些是真的前提下，进一步提高压缩比，以换取更高的编码效率，这就需要一些标准来对图像压缩的质量进行评价。在我们的程序中，主要运用的是均方根误差来对图像压缩的质量进行评价。

 

四、参考文献及资料

【1】ShoichiroNakamura .科学计算引论.电子工业出版社，2002.

【2】低通滤波技术在数字图像处理中的应用以及仿真实现（百度资料）.

【3】语音信号的滤波处理（课程设计）.

【4】胡乡峰，卫金茂. 基于奇异值分解的图像压缩.东北师范大学学报，2006.

【5】曾超，张卫东.基于奇异值分解的图像压缩及其Matlab实现.高校理科研究，2010.