一个用MATLAB自带函数SVD做低秩压缩和去噪的小例子

首先放一下奇异值分解的原理：
奇异值分解的揭秘（一）：矩阵的奇异值分解过程

下面是用用MATLAB自带函数SVD做低秩压缩和去噪的实际过程：

1.打开原图

（以前一直认为clc和clear都是在调试的时候在命令行里用，最近发现了在编辑器里的必要之处，保证多次调试的时候不会出现，如果注释掉某段代码而自己忘记了，就不会沿用之前的变量值出现匪夷所思的场面但是自己发现不了实际情况，而一再认为参数有问题的情况。）

clc;
clear;

img_name = 'freckle1.jpg';%原图：带有雀斑的人脸
img = imread(img_name);

（imread是把图片读进来，如果不读就会报错，是不能直接打开的，imshow是展示图片）

gray_img = rgb2gray(img);%将图片改为灰度图
img_gray = double(gray_img);%将图片改为双精度型用于计算
[U,S,V] = svd(img_gray);

必须转换为double类型才能计算

注意svd不能处理rgb图像，直接用会报错
（奇异值分解只能针对二维矩阵或高维矩阵降为二维后使用）

>> [U,S,V] = svd(orig_img);
错误使用 svd
输入必须为 2 维。

2.计算要压缩的比例

%计算对角矩阵S的非零元素个数
S_diag = diag(S);
% A=find(S_diag);
% size(A);
k = sum(S_diag~=0);
k1=fix(0.05*k);k2=fix(0.1*k);k3=fix(0.3*k);%分别输出压缩为0.05，0.1和0.3时候的图像

在svd中矩阵S为对角阵，其非零元素个数即为矩阵的秩（奇异值个数），通过sum函数可以直接求出秩，或者用find函数找出不等于0的项形成数组再求size

3.显示

要先转换为uint8类型才能显示

%转换回uint类型用于显示
img_svd1 = uint8(BD1);
img_svd2 = uint8(BD2);
img_svd3 = uint8(BD3);

画图

%画图
subplot(2,2,1);
imshow(gray_img);
title('原图');
subplot(2,2,2);
imshow(img_svd1);
title(['压缩至0.05,相似度为',num2str(BD1_error_mean)]);
subplot(2,2,3);
imshow(img_svd2);
title(['压缩至0.1,相似度为',num2str(BD2_error_mean)]);
subplot(2,2,4);
imshow(img_svd3);
title(['压缩至0.3,相似度为',num2str(BD3_error_mean)]);

求与原图的近似度过程如下

4.求与原图的近似度

将压缩后的U,S,V矩阵乘起来得到压缩后的矩阵，将原图每一个像素点减去压缩后像素点比上原图像素点的值求绝对值，再求平均数得到error，1-error近似为与原图的相似度

%用1-error近似表示相似度
BD1_error=mean(abs((img_gray-BD1)./img_gray));%取前5%个特征值时的error
BD1_error2=find(BD1_error==Inf);%如果原图灰度值为0，就会出现Inf影响求均值，所以将这些Inf值变为0
BD1_error(BD1_error2)=0;
BD1_error_mean=1-mean(BD1_error);

注意如果原图中有0就会出现Inf影响计算，所以可以除去Inf值再做平均（不可以先除去0值因为会有较大误差）
如果要求多个压缩比的近似度其它可以不必重新寻找Inf值因为原图相同但必须分别将Inf值置0

运行结果

可以看出，取5%的奇异值时近似度已经达到97%以上，结果非常好

6.降噪

img=double(gray_img) +wgn(size(gray_img,1),size(gray_img,2),20);%加高斯白噪声
img_show=uint8(img);%用来展示的噪声图

这里选用了一种MATLAB自带的加高斯白噪声的方法，以下链接还有很多非常好的方法，这里就不展开描述了
matlab给数字图像加高斯白噪声的几种方法

运行结果

这里只加了20dB的白噪声，但是可以看出效果不是特别理想，并不能非常接近原图如果噪声增强效果会更差，所以笔者认为svd在低秩压缩表现比降噪要好很多