% 基于SVD的图像压缩演示
clear; clc; close all;
%% 1. 读取并显示原始图像
try
% 尝试读取内置图像
original_img = imread('ex2.jpg');
catch
% 如果内置图像不存在，创建一个示例图像
fprintf('使用内置图像失败，创建示例图像...\n');
[X, Y] = meshgrid(1:256, 1:256);
original_img = uint8(128 + 100 * sin(X/20) .* cos(Y/15));
end
% 如果图像是彩色的，转换为灰度图
if ndims(original_img) == 3 && size(original_img, 3) == 3
original_img = rgb2gray(original_img);
end
% 转换为double类型以便计算
A = double(original_img);
% 显示原始图像
figure('Position', [100, 100, 1200, 600]);
subplot(2, 3, 1);
imshow(original_img);
title('原始图像', 'FontSize', 12);
xlabel(['尺寸: ', num2str(size(A, 1)), '×', num2str(size(A, 2))]);
%% 2. 对图像矩阵进行SVD分解
fprintf('正在进行SVD分解...\n');
[U, S, V] = svd(A);
% 获取奇异值
singular_values = diag(S);
total_singular_values = length(singular_values);
fprintf('SVD分解完成！总奇异值个数: %d\n', total_singular_values);
%% 3. 使用不同数量的奇异值重建图像
k_values = [1, 5, 20, 50, 100]; % 选择不同的k值进行测试
% 计算原始图像需要存储的元素数量
original_elements = numel(A);
fprintf('原始图像需要存储的元素数量: %d\n', original_elements);
for i = 1:length(k_values)
k = k_values(i);
if k > total_singular_values
k = total_singular_values;
end
% 使用前k个奇异值重建图像
A_compressed = U(:, 1:k) * S(1:k, 1:k) * V(:, 1:k)';
% 确保数值在合理范围内
A_compressed = max(min(A_compressed, 255), 0);
% 计算压缩比
compressed_elements = size(U, 1)*k + k + size(V, 1)*k;
compression_ratio = compressed_elements / original_elements;
% 显示重建后的图像
subplot(2, 3, i+1);
imshow(uint8(A_compressed));
title_str = sprintf('k = %d\n压缩比: %.1f%%', k, compression_ratio*100);
title(title_str, 'FontSize', 10);
% 计算PSNR（峰值信噪比）
mse = mean((A(:) - A_compressed(:)).^2);
if mse > 0
psnr_val = 20*log10(255/sqrt(mse));
else
psnr_val = Inf;
end
xlabel(sprintf('PSNR: %.1f dB', psnr_val));
end
%% 4. 绘制奇异值分析图
figure('Position', [100, 500, 1000, 400]);
% 绘制奇异值衰减曲线
subplot(1, 2, 1);
plot_range = min(100, total_singular_values);
plot(1:plot_range, singular_values(1:plot_range), 'b-', 'LineWidth', 2);
xlabel('奇异值序号');
ylabel('奇异值大小');
title('前100个奇异值');
grid on;
% 绘制累积能量
subplot(1, 2, 2);
cumulative_energy = cumsum(singular_values.^2) / sum(singular_values.^2);
plot(1:plot_range, cumulative_energy(1:plot_range), 'r-', 'LineWidth', 2);
xlabel('使用的奇异值数量 (k)');
ylabel('累积能量比例');
title('信息保留程度 vs k值');
grid on;
% 标记关键点
hold on;
marker_k = [1, 5, 20, 50];
for j = 1:length(marker_k)
k = marker_k(j);
if k <= plot_range
plot(k, cumulative_energy(k), 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'red');
text(k, cumulative_energy(k)-0.1, sprintf('k=%d\n%.1f%%', k, cumulative_energy(k)*100), ...
'HorizontalAlignment', 'center', 'FontSize', 8);
end
end
%% 5. 交互式压缩演示
figure('Position', [200, 200, 800, 600]);
fprintf('\n=== 交互式压缩演示 ===\n');
fprintf('请输入想要的压缩级别 (k值，1-%d)，输入0退出: ', total_singular_values);
while true
k_input = input('');
if k_input == 0
break;
end
if k_input < 1 || k_input > total_singular_values
fprintf('请输入1-%d之间的数字: ', total_singular_values);
continue;
end
% 使用指定k值重建图像
A_custom = U(:, 1:k_input) * S(1:k_input, 1:k_input) * V(:, 1:k_input)';
A_custom = max(min(A_custom, 255), 0);
% 计算指标
compressed_elems = size(U, 1)*k_input + k_input + size(V, 1)*k_input;
comp_ratio = compressed_elems / original_elements;
info_kept = cumulative_energy(k_input) * 100;
mse_custom = mean((A(:) - A_custom(:)).^2);
if mse_custom > 0
psnr_custom = 20*log10(255/sqrt(mse_custom));
else
psnr_custom = Inf;
end
% 显示结果
clf;
subplot(1, 2, 1);
imshow(original_img);
title('原始图像', 'FontSize', 14);
subplot(1, 2, 2);
imshow(uint8(A_custom));
title_str = sprintf('压缩图像 (k=%d)', k_input);
title(title_str, 'FontSize', 14);
% 显示压缩信息
info_str = sprintf(['压缩比: %.2f%%\n', ...
'信息保留: %.1f%%\n', ...
'PSNR: %.1f dB\n', ...
'存储元素: %d → %d'], ...
comp_ratio*100, info_kept, psnr_custom, ...
original_elements, compressed_elems);
text(10, 30, info_str, 'Color', 'white', 'FontSize', 12, ...
'BackgroundColor', 'black', 'VerticalAlignment', 'top');
fprintf('k=%d: 压缩比=%.2f%%, 信息保留=%.1f%%, PSNR=%.1f dB\n', ...
k_input, comp_ratio*100, info_kept, psnr_custom);
fprintf('继续输入k值(1-%d)或输入0退出: ', total_singular_values);
end
%% 6. 压缩效果总结
fprintf('\n=== 压缩效果总结 ===\n');
fprintf('k值\t压缩比\t\t信息保留\tPSNR(dB)\n');
fprintf('------------------------------------------------\n');
for k = [1, 5, 10, 20, 50, 100]
if k <= total_singular_values
compressed_elems = size(U, 1)*k + k + size(V, 1)*k;
comp_ratio = compressed_elems / original_elements;
info_kept = cumulative_energy(k) * 100;
A_temp = U(:, 1:k) * S(1:k, 1:k) * V(:, 1:k)';
mse_temp = mean((A(:) - A_temp(:)).^2);
if mse_temp > 0
psnr_temp = 20*log10(255/sqrt(mse_temp));
else
psnr_temp = Inf;
end
fprintf('%d\t%.2f%%\t\t%.1f%%\t\t%.1f\n', k, comp_ratio*100, info_kept, psnr_temp);
end
end
fprintf('\n程序运行完成！\n');
公众号