使用SVD处理图像模拟演示

参考文献

https://github.com/matzewolf/Image_compression_SVD/blob/master/svd_compress.m

MATLAB代码

clc;
clearvars;
close all;A_org=imread("lena256.bmp");
compr=20;
A_org=double(A_org);A_red  = svd_compress( A_org, compr );
subplot(1,2,1),imshow(A_org,[]);
subplot(1,2,2),imshow(A_red,[]);function [ A_red ] = svd_compress( A_org, compr )% svd_compress compresses an input matrix (e.g. an image) using the
% Singular Value Decomposition (SVD).
%   Input args: A_org: Any matrix with double real entries, e.g. an image 
%   file (converted from uint8 to double).
%   compr: Quality of compression. If 0 <= compr < 1, it only keeps
%   Singular Values (SVs) larger than compr times the biggest SV. If 1 <= 
%   compr <= number of SVs, it keeps the biggest compr SVs. Otherwise the 
%   function returns an error.
%   Output args: A_red: Compressed version of A_org in double using the
%   SVD, e.g. an image file (convert from double to uint8).% SVD on the original matrix
[U,S,V] = svd(A_org);% Extract Singular Values (SVs)
singvals = diag(S);% Determine SVs to be saved
if compr >= 0 && compr < 1% only SVs bigger than compr times biggest SVindices = find(singvals >= compr * singvals(1));
elseif compr >= 1 && compr <= length(singvals)% only the biggest compr SVsindices = 1:compr;
else% return errorerror('Incorrect input arg: compr must satisfy 0 <= compr <= number of Singular Values');
end% Truncate U,S,V
U_red = U(:,indices);
S_red = S(indices,indices);
V_red = V(:,indices);% Calculate compressed matrix
A_red = U_red * S_red * V_red';end

运行结果

% Image Compression with Singular Value Decomposition (SVD).
%   This script uses the SVD for Image Compression, analyses the algorithm
%   (also with Information Theory) and visualizes the results.close all; clear; clc;
tic;
COL = 256; % number of colors in uint8, so 2^8 = 256.%% Compression% Original image matrix
Lena_org = imread('Lena.bmp'); % in uint8
Lena = double(Lena_org); % in double% Call compressing function (and measure performance)
compr = 0.01; % change compr to change quality
tic;
Lena_red = uint8(svd_compress(Lena,compr));
func_time = toc; % compression function execution time
fprintf('Execution time of svd_compress: %d seconds.\n',func_time);% Save compressed image
imwrite(Lena_red,'ReducedLena.bmp');%% Analysis of the algorithm% SVD on the image
[U,S,V] = svd(Lena);% Extract Singular Values (SVs)
singvals = diag(S);% Determine SVs to be saved
if compr >= 0 && compr < 1% only SVs bigger than compr times biggest SVindices = find(singvals >= compr * singvals(1));
elseif compr >= 1 && compr <= length(singvals)% only the biggest compr SVsindices = 1:compr;
else% return errorerror(...'Incorrect input arg: compr must satisfy 0 <= compr <= number of Singular Values');
end% Size of the image
m = size(Lena,1);
n = size(Lena,2);
storage = m*n;
fprintf('Size of image: %d px by %d px, i.e. uses %d px of storage.\n',m,n,storage);% SVs and reduced storage
r = min([m,n]); % original number of SVs
r_red = length(indices); % to be saved number of SVs
r_max = floor(m*n/(m+n+1)); % maximum to be saved number of SVs for compression
storage_red = m*r_red + n*r_red + r_red;
if compr >= 0 && compr < 1% only SVs bigger than compr times biggest SVfprintf('The smallest SV chosen to be smaller than %d of the biggest SV.\n',compr);
elseif compr >= 1 && compr <= length(singvals)% only the biggest compr SVs
else% return errorfprintf('There was some error before. Analysis cannot continue.\n')
end
fprintf('Out of %d SVs, only %d SVs saved ',r,r_red);
fprintf('(Maximum number of SVs for compression: %d SVs).\n',r_max);
fprintf('Reduced storage: %d px.\n',storage_red);% Determine made error
error = 1 - sum(singvals(indices))/sum(singvals);
fprintf('Made error: %d.\n',error);
errorImage = Lena_org - Lena_red;% Entropy
entropy_org = entropy(Lena_org);
fprintf('Entropy of original image: %d bit.\n',entropy_org);
entropy_red = entropy(Lena_red);
fprintf('Entropy of compressed image: %d bit.\n',entropy_red);
entropy_err = entropy(errorImage);
fprintf('Entropy of error image: %d bit.\n',entropy_err);% 1D Histogram: Original Probability
[orgProb,~,~] = histcounts(Lena_org,1:(COL+1),'Normalization','probability');% 2D Histogram: Joint Probabiltiy
[jointProb,~,~] = histcounts2(Lena_red,Lena_org,...1:(COL+1),1:(COL+1),'Normalization','probability');% Joint Entropy
p_logp_nan = jointProb.*log2(jointProb);
p_logp = p_logp_nan(isfinite(p_logp_nan));
joint_entropy = -sum(p_logp);
fprintf('Joint entropy: %d bit.\n',joint_entropy);% Mutual Information
mi = entropy_org + entropy_red - joint_entropy;
fprintf('Mutual information: %d bit.\n',mi);% Conditional Probability
condProb = jointProb./orgProb;
condProb(isnan(condProb)|isinf(condProb))=0; % all NaN and inf converted to zero
col_sum = sum(condProb,1); % test if condProb really sums up to 1 columnwise%% Relationship between selcted SVs and ...numSVals = 1:1:r; %SVs for which the properties are calculated% ...used storage
storageSV = m*numSVals + n*numSVals + numSVals;% ...made error and entropies (compressed and error)
displayedError = zeros(size(numSVals));
entropySV = zeros(4,length(numSVals));% 1st row entropy of compressed image, 2nd row entropy of error image% 3rd row joint entropy, 4th row mutual information
j = 1; % position in the display vectors
for i = numSVals% store S in a temporary matrixS_loop = S;% truncate SS_loop(i+1:end,:) = 0;S_loop(:,i+1:end) = 0;% construct Image using truncated SLena_red_loop = uint8(U*S_loop*V');% construct error imageLena_err_loop = Lena_org - Lena_red_loop;% compute errorerror_loop = 1 - sum(diag(S_loop))/sum(diag(S));% add error to display vectordisplayedError(j) = error_loop;% compute entropy of compressed image and add to row 1 of display matrixentropySV(1,j) = entropy(Lena_red_loop);% compute entropy of error image and add to row 2 of display matrixentropySV(2,j) = entropy(Lena_err_loop);% compute joint entropy of original and compresed image[jointProb_loop,~,~] = histcounts2(Lena_org,Lena_red_loop,[COL COL],...'Normalization','probability');p_logp_nan_loop = jointProb_loop.*log2(jointProb_loop);p_logp_loop = p_logp_nan_loop(isfinite(p_logp_nan_loop));entropySV(3,j) = -sum(p_logp_loop);% compute mutual information of original and compressed imageentropySV(4,j) = entropy_org + entropySV(1,j) - entropySV(3,j);% update positionj = j + 1;
end%% Figure 1fig1 = figure('Name','Images and Histograms',...'units','normalized','outerposition',[0 0 1 1]);% Original image
subplot(2,3,1)
imshow(uint8(Lena))
title('Original image')% Histogram of original image
subplot(2,3,4)
imhist(Lena_org)
title('Histogram of original image')% Compressed image
subplot(2,3,2)
imshow(uint8(Lena_red))
title('Compressed image')% Histogram of compressed image
subplot(2,3,5)
imhist(Lena_red)
title('Histogram of compressed image')% Error image
subplot(2,3,3)
imshow(uint8(errorImage))
title('Error image')% Histogram of error image
subplot(2,3,6)
imhist(errorImage)
title('Histogram of error image')%% Figure 2fig2 = figure('Name','Joint Histogram',...'units','normalized','outerposition',[0 0 1 1]);% 2D Histogram: Joint PDF
histogram2(Lena_red,Lena_org,1:(COL+1),1:(COL+1),...'Normalization','probability','FaceColor','flat')
colorbar
title('Joint Histogram')
xlabel('Compressed image')
ylabel('Original image')
zlabel('Joint Probability')%% Figure 3fig3 = figure('Name','Properties over selected Singular Values',...'units','normalized','outerposition',[0 0 1 1]);% Used storage over saved SVs
subplot(2,2,1)
plot(numSVals, storage.*ones(size(numSVals))) % original storage (horizontal)
hold on
plot(numSVals, storageSV)
legend('Original storage', 'Storage of SVD','Location','northwest')
xlabel('Number of saved Singular Values')
ylabel('Used storage [px]')
title('Used storage over saved SVs')% Compression error over saved SVs
subplot(2,2,3)
plot(numSVals, displayedError)
xlabel('Number of saved Singular Values')
ylabel('Compression error [-]')
title('Compression error over saved SVs')% Entropies over saved SVs
subplot(2,2,[2,4])
plot(numSVals, entropy_org.*ones(size(numSVals))) % original entropy (horizontal)
hold on
plot(numSVals, entropySV)
legend('Original entropy', 'Compression entropy', 'Error entropy',...'Joint entropy','Mutual information','Location','southoutside')
xlabel('Number of saved Singular Values')
ylabel('Entropies [bit]')
title('Entropies over saved SVs')%% Save figuressaveas(fig1, 'Results.png');
saveas(fig2, 'Joint_Histogram.png');
saveas(fig3, 'Analysis.png');%% Execution timeexecution_time = toc; % total script execution time
fprintf('Total execution time of svd_lena_script: %d seconds.\n',execution_time);

运行结果