使用SVD在MATLAB中压缩图像

Question

我是一个Matlab新手，但正在尝试编写一些针对灰度图像的图像压缩代码。

问题

如何使用SVD修剪低值特征值以重构压缩图像？

目前的工作/尝试

到目前为止，我的代码是：

B=imread('images1.jpeg');   
B=rgb2gray(B);  
doubleB=double(B);  
%read the image and store it as matrix B, convert the image to a grayscale
photo and convert the matrix to a class 'double' for values 0-255  
[U,S,V]=svd(doubleB);

这使我能够成功地使用存储在变量S中的特征值分解图像矩阵。

如何截断S（大小为167x301，类别为double）？假设我只想取其中前100个（或任何n），我该如何做并重建压缩的图像？

更新的代码/想法

与其在评论部分放置一堆代码，这是我目前拥有的草稿。我已经成功创建了压缩图像，手动更改N的值，但我还想做两件事：

1- 显示各种压缩的图像面板（比如，对于N = 5、10、25等运行一个循环）。

2- 以某种方式计算每个图像与原始图像之间的差异（误差），并将其绘制成图形。

我不擅长理解循环和输出，但这是我尝试过的内容：

B=imread('images1.jpeg');  
B=rgb2gray(B);  
doubleB=im2double(B);%  
%read the image and store it as matrix B, convert the image to a grayscale  
%photo and convert the image to a class 'double'  
[U,S,V]=svd(doubleB);   
C=S;  
for N=[5,10,25,50,100]  
C(N+1:end,:)=0;  
C(:,N+1:end)=0;  
D=U*C*V';  
%Use singular value decomposition on the image doubleB, create a new matrix  
%C (for Compression diagonal) and zero out all entries above N, (which in  
%this case is 100). Then construct a new image, D, by using the new  
%diagonal matrix C.  
imshow(D);  
error=C-D;  
end

显然存在一些错误，因为我没有得到多张图片，也不知道如何“绘制”误差矩阵。

- Justin

---

更多有关DSP的信息请参见：https://dsp.stackexchange.com/questions/18673/image-compression-using-the-svd-in-matlab。 - GeorgeIrwin

Answer 1

尽管这个问题有点老了，但它帮助我很多地理解了SVD。我修改了你在问题中写的代码以使其能够工作。

我相信你可能已经解决了这个问题，但是为了方便日后访问本页面的任何人，我在此处包含了完整的代码以及输出图像和图表。

以下是代码：

close all
clear all
clc

%reading and converting the image
inImage=imread('fruits.jpg');
inImage=rgb2gray(inImage);
inImageD=double(inImage);

% decomposing the image using singular value decomposition
[U,S,V]=svd(inImageD);

% Using different number of singular values (diagonal of S) to compress and
% reconstruct the image
dispEr = [];
numSVals = [];
for N=5:25:300
    % store the singular values in a temporary var
    C = S;

    % discard the diagonal values not required for compression
    C(N+1:end,:)=0;
    C(:,N+1:end)=0;

    % Construct an Image using the selected singular values
    D=U*C*V';


    % display and compute error
    figure;
    buffer = sprintf('Image output using %d singular values', N)
    imshow(uint8(D));
    title(buffer);
    error=sum(sum((inImageD-D).^2));

    % store vals for display
    dispEr = [dispEr; error];
    numSVals = [numSVals; N];
end

% dislay the error graph
figure; 
title('Error in compression');
plot(numSVals, dispEr);
grid on
xlabel('Number of Singular Values used');
ylabel('Error between compress and original image');

将以下内容应用于图像： 只使用前5个奇异值得到以下结果: 只使用前30个奇异值得到以下结果: 只使用前55个奇异值得到以下结果: 随着奇异值数量的增加，误差变化如下图所示： 可以看出，使用大约前200个奇异值可以使误差近似为零。

- masad

Answer 2

首先，我假设您已经意识到SVD并不是单个图像中去相关像素的最佳工具。 但这是一个很好的练习。

好的，我们知道B = U*S*V'。 我们知道S是对角线矩阵，并按大小排序。 因此，仅使用前几个S的值，您将获得您的图像的近似值。 假设C=U*S2*V'，其中S2是您修改后的S。 U和V的大小未更改，因此现在最简单的方法就是将您不想使用的S元素归零并运行重建。(这样做的最简单方法是：S2=S; S2(N+1:end, :) = 0; S2(:, N+1:end) = 0;)。

现在是压缩部分。U和V都是完整的，因此无论S2发生什么变化，数据量都不会改变。 但是看看U*S2会发生什么。（绘制图像）。 如果您保留了S2中的N个奇异值，则S2的前N行是唯一的非零行。 压缩！ 但是您仍然必须处理V。 一旦您完成了(U*S2)，就无法再使用相同的技巧，因为比S2本身更多的U*S2元素是非零的。 我们如何在两侧都使用S2？ 好吧，它是对角线矩阵，所以使用D=sqrt(S2)，现在C=U*D*D*V'。 因此，U*D只有N个非零行，D*V'只有N个非零列。 仅传输那些量，您就可以重建近似于B的C。

- Peter

---

非常感谢您的详细解释。我会仔细查看并在有问题或困难时再回来咨询。谢谢！ - Justin

目前我正在使用以下代码：B=imread('images1.jpeg'); B=rgb2gray(B); doubleB=double(B); [U,S,V]=svd(doubleB); C=S; N=100; C(N+1:end,:)=0; C(:,N+1:end)=0; D=U*C*V'; imshow(D); 无论我将N设置为多少，新图像似乎都是相同的（并且看起来非常草图）。 参考一下，S是167x301。 我做错了什么？ - Justin

嗯，实际上，当我运行imshow(doubleB)时，图像看起来非常糟糕...与原始图像完全不同。 - Justin

尝试使用im2double(B)代替double(B)。 - Peter

我刚刚在原始帖子中添加了一些额外的代码和想法。我遇到了一些困难，不知道如何从循环中输出图像（改变N），并绘制压缩和原始之间的差异/误差。 - Justin

请查看subplot以显示多个图像。尝试使用imagesc而不是image来显示错误图像。并且/或者将像素误差图像平方并取平均值，以获得每个重建的单个误差指标（均方误差）。将该指标存储在另一个向量中，以便您可以绘制N vs. MSE。 - Peter

Answer 3

例如，这是一个512 x 512的黑白图像 Lena:

我们计算 Lena 的 SVD。选择最大奇异值的1%以上的奇异值，我们只剩下53个奇异值。使用这些奇异值和相应的（左右）奇异向量重构 Lena，我们得到了 Lena 的低秩近似：

我们可以存储2 x (512 x 53) + 53 = 54325个值，而不是存储512 * 512 = 262144个值（每个值占8位），这大约是原始大小的20%。这是SVD用于进行有损图像压缩的示例之一。

---

以下是MATLAB代码：

% open Lena image and convert from uint8 to double
Lena = double(imread('LenaBW.bmp'));

% perform SVD on Lena
[U,S,V] = svd(Lena);

% extract singular values
singvals = diag(S);

% find out where to truncate the U, S, V matrices
indices = find(singvals >= 0.01 * singvals(1));

% reduce SVD matrices
U_red = U(:,indices);
S_red = S(indices,indices);
V_red = V(:,indices);

% construct low-rank approximation of Lena
Lena_red = U_red * S_red * V_red';

% print results to command window
r = num2str(length(indices));
m = num2str(length(singvals));
disp(['Low-rank approximation used ',r,' of ',m,' singular values']);

% save reduced Lena
imwrite(uint8(Lena_red),'Reduced Lena.bmp');

- Rodrigo de Azevedo

Answer 4

取前n个最大的特征值及其对应的特征向量可能会解决您的问题。对于PCA，原始数据乘以前n个升序特征向量将构建您的图像，其中n x d表示特征向量的数量。

- berkay