% 基于SVD的图像压缩演示

clear; clc; close all;

%% 1. 读取并显示原始图像

% 读取图像（请将 'your_image.jpg' 替换为你的图像路径）

original_img = imread('ex2.jpg'); 

% 如果图像是彩色的，转换为灰度图

if size(original_img, 3) == 3

    original_img = rgb2gray(original_img);

end

% 转换为double类型以便计算

A = double(original_img);

% 显示原始图像

figure('Position', [100, 100, 1200, 800]);

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);

    % 使用前k个奇异值重建图像

    A_compressed = U(:, 1:k) * S(1:k, 1:k) * V(:, 1:k)';

    % 计算压缩比

    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(sprintf('k = %d (压缩比: %.1f%%)', k, compression_ratio*100), 'FontSize', 10);

    xlabel(sprintf('需要存储: %d个元素', compressed_elements));

    % 计算并显示PSNR（峰值信噪比）来衡量重建质量

    mse = mean((A(:) - A_compressed(:)).^2);

    psnr_val = 20*log10(255/sqrt(mse));

    text(10, 20, sprintf('PSNR: %.1f dB', psnr_val), 'Color', 'white', 'FontSize', 8, 'BackgroundColor', 'black');

end

%% 4. 绘制奇异值衰减曲线

figure('Position', [100, 500, 800, 400]);

% 绘制奇异值

subplot(1, 2, 1);

semilogy(1:min(100, total_singular_values), singular_values(1:min(100, total_singular_values)), '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:min(100, total_singular_values), cumulative_energy(1:min(100, total_singular_values)), 'r-', 'LineWidth', 2);

xlabel('使用的奇异值数量 (k)');

ylabel('累积能量比例');

title('信息保留程度 vs k值');

grid on;

% 标记几个关键点

hold on;

marker_k = [1, 5, 20, 50];

for k = marker_k

    if k <= length(cumulative_energy)

        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. 显示压缩效果对比

fprintf('\n=== 压缩效果总结 ===\n');

fprintf('k值\t压缩比\t\t信息保留\n');

for k = [1, 5, 20, 50, 100]

    if k <= total_singular_values

        compressed_elements = size(U, 1)*k + k + size(V, 1)*k;

        compression_ratio = compressed_elements / original_elements;

        info_retained = cumulative_energy(k) * 100;

        fprintf('%d\t%.2f%%\t\t%.1f%%\n', k, compression_ratio*100, info_retained);

    end

end

fprintf('\n程序运行完成！\n');

% 基于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');