基于低秩约束去除图像中稀疏噪声是一个经典的图像处理问题，通常使用鲁棒主成分分析（Robust PCA, RPCA）方法。
1. 理论基础
问题建模
图像去噪问题可以表示为：
M = L + S + N
其中：
- M：观测到的含噪图像
- L：低秩成分（干净图像）
- S：稀疏噪声（如椒盐噪声、遮挡等）
- N：高斯噪声
数学模型
RPCA优化问题：
其中：
‖L‖_*：核范数（奇异值之和），促进低秩‖S‖_1：L1范数，促进稀疏性λ：权衡参数
参考代码 使用图像的低秩约束去除图像中稀疏噪声 www.youwenfan.com/contentcsl/81497.html
2. MATLAB实现方法
方法一：使用主成分追踪（Principal Component Pursuit）
function [L, S] = rpca_denoise(M, lambda, mu, tol, max_iter)
% RPCA图像去噪
% 输入:
% M - 含噪图像矩阵
% lambda - 稀疏项权重 (默认: 1/sqrt(max(size(M))))
% mu - 增强拉格朗日参数 (默认: 1.25/norm(M,2))
% tol - 收敛容差 (默认: 1e-7)
% max_iter - 最大迭代次数 (默认: 1000)
% 输出:
% L - 低秩成分 (去噪后图像)
% S - 稀疏噪声
if nargin < 2 || isempty(lambda)
lambda = 1 / sqrt(max(size(M)));
end
if nargin < 3 || isempty(mu)
mu = 1.25 / norm(M, 2);
end
if nargin < 4 || isempty(tol)
tol = 1e-7;
end
if nargin < 5 || isempty(max_iter)
max_iter = 1000;
end
[m, n] = size(M);
% 初始化变量
L = zeros(m, n);
S = zeros(m, n);
Y = zeros(m, n);
for iter = 1:max_iter
% 更新L (低秩成分)
[U, sigma, V] = svd(M - S + Y/mu, 'econ');
sigma = diag(sigma);
svp = length(find(sigma > 1/mu));
if svp >= 1
sigma = sigma(1:svp) - 1/mu;
else
svp = 1;
sigma = 0;
end
L = U(:, 1:svp) * diag(sigma) * V(:, 1:svp)';
% 更新S (稀疏噪声)
temp = M - L + Y/mu;
S = sign(temp) .* max(abs(temp) - lambda/mu, 0);
% 更新拉格朗日乘子
Z = M - L - S;
Y = Y + mu * Z;
% 检查收敛性
err = norm(Z, 'fro') / norm(M, 'fro');
if mod(iter, 50) == 0
fprintf('迭代 %d, 误差: %.2e\n', iter, err);
end
if err < tol
fprintf('在迭代 %d 收敛\n', iter);
break;
end
end
if iter == max_iter
fprintf('达到最大迭代次数\n');
end
end
方法二：加速近端梯度方法
function [L, S] = apg_rpca(M, lambda, max_iter, tol)
% 加速近端梯度方法求解RPCA
if nargin < 2
lambda = 1 / sqrt(max(size(M)));
end
if nargin < 3
max_iter = 500;
end
if nargin < 4
tol = 1e-6;
end
[m, n] = size(M);
% 初始化
L = zeros(m, n);
L_prev = L;
S = zeros(m, n);
t = 1;
t_prev = 1;
for iter = 1:max_iter
% 计算加速点
Y_L = L + ((t_prev - 1)/t) * (L - L_prev);
% 梯度步
G = Y_L - 0.5 * (Y_L + S - M);
% 更新L (奇异值阈值)
[U, sigma, V] = svd(G, 'econ');
sigma = diag(sigma);
tau = 0.5;
svp = length(find(sigma > tau));
if svp >= 1
sigma = sigma(1:svp) - tau;
else
svp = 1;
sigma = 0;
end
L_prev = L;
L = U(:, 1:svp) * diag(sigma) * V(:, 1:svp)';
% 更新S (软阈值)
Y_S = S + ((t_prev - 1)/t) * (S - S);
G_S = Y_S - 0.5 * (L + Y_S - M);
S = sign(G_S) .* max(abs(G_S) - lambda * 0.5, 0);
% 更新加速参数
t_prev = t;
t = (1 + sqrt(1 + 4 * t^2)) / 2;
% 检查收敛性
residual = norm(M - L - S, 'fro') / norm(M, 'fro');
if mod(iter, 50) == 0
fprintf('迭代 %d, 残差: %.2e\n', iter, residual);
end
if residual < tol
fprintf('收敛于迭代 %d\n', iter);
break;
end
end
end
3. 完整的图像去噪示例
function low_rank_denoising_demo()
% 低秩约束图像去噪完整示例
%% 1. 读取和准备数据
% 读取图像
clean_img = imread('lena.png');
if size(clean_img, 3) > 1
clean_img = rgb2gray(clean_img);
end
clean_img = im2double(clean_img);
% 添加混合噪声
noisy_img = add_mixed_noise(clean_img, 0.1, 0.05);
%% 2. 参数设置
[m, n] = size(noisy_img);
lambda = 1 / sqrt(max(m, n));
max_iter = 300;
tol = 1e-6;
%% 3. 使用RPCA去噪
fprintf('开始RPCA去噪...\n');
tic;
[L_denoised, S_noise] = rpca_denoise(noisy_img, lambda, [], tol, max_iter);
time_elapsed = toc;
fprintf('去噪完成，耗时: %.2f 秒\n', time_elapsed);
%% 4. 结果评估
% 计算性能指标
psnr_original = psnr(noisy_img, clean_img);
psnr_denoised = psnr(L_denoised, clean_img);
ssim_original = ssim(noisy_img, clean_img);
ssim_denoised = ssim(L_denoised, clean_img);
fprintf('性能指标:\n');
fprintf('含噪图像 PSNR: %.2f dB, SSIM: %.4f\n', psnr_original, ssim_original);
fprintf('去噪图像 PSNR: %.2f dB, SSIM: %.4f\n', psnr_denoised, ssim_denoised);
%% 5. 结果显示
figure('Position', [100, 100, 1200, 800]);
subplot(2, 3, 1);
imshow(clean_img); title('原始干净图像');
subplot(2, 3, 2);
imshow(noisy_img); title('含噪图像');
subplot(2, 3, 3);
imshow(L_denoised); title('RPCA去噪结果');
subplot(2, 3, 4);
imshow(S_noise, []); title('提取的稀疏噪声');
subplot(2, 3, 5);
% 显示残差
residual = abs(clean_img - L_denoised);
imshow(residual, []); title('去噪残差');
colorbar;
subplot(2, 3, 6);
% 奇异值分布比较
sv_clean = svd(clean_img);
sv_noisy = svd(noisy_img);
sv_denoised = svd(L_denoised);
semilogy(sv_clean(1:50), 'b-', 'LineWidth', 2); hold on;
semilogy(sv_noisy(1:50), 'r--', 'LineWidth', 2);
semilogy(sv_denoised(1:50), 'g-.', 'LineWidth', 2);
legend('干净图像', '含噪图像', '去噪图像');
title('奇异值分布');
xlabel('索引'); ylabel('奇异值');
grid on;
end
function noisy_img = add_mixed_noise(clean_img, salt_pepper_ratio, gaussian_var)
% 添加混合噪声
noisy_img = imnoise(clean_img, 'salt & pepper', salt_pepper_ratio);
noisy_img = imnoise(noisy_img, 'gaussian', 0, gaussian_var);
end
4. 针对彩色图像的处理
function [L_color, S_color] = rpca_color_denoise(color_img, lambda)
% 彩色图像RPCA去噪
[h, w, c] = size(color_img);
% 将彩色图像重塑为矩阵 (每个通道作为一列)
M_matrix = reshape(color_img, h*w, c);
% 对彩色通道应用RPCA
[L_matrix, S_matrix] = rpca_denoise(M_matrix, lambda);
% 重塑回图像格式
L_color = reshape(L_matrix, h, w, c);
S_color = reshape(S_matrix, h, w, c);
% 确保值在[0,1]范围内
L_color = max(0, min(1, L_color));
S_color = max(-1, min(1, S_color));
end
5. 基于块的低秩去噪
function denoised_img = patch_based_low_rank_denoise(noisy_img, patch_size, overlap)
% 基于图像块的局部低秩去噪
if nargin < 2
patch_size = 16;
end
if nargin < 3
overlap = 4;
end
[h, w] = size(noisy_img);
denoised_img = zeros(h, w);
weight = zeros(h, w);
step = patch_size - overlap;
for i = 1:step:h-patch_size+1
for j = 1:step:w-patch_size+1
% 提取图像块
patch = noisy_img(i:i+patch_size-1, j:j+patch_size-1);
% 对块应用RPCA
lambda = 1 / sqrt(patch_size);
[L_patch, ~] = rpca_denoise(patch, lambda);
% 累积到输出图像
denoised_img(i:i+patch_size-1, j:j+patch_size-1) = ...
denoised_img(i:i+patch_size-1, j:j+patch_size-1) + L_patch;
weight(i:i+patch_size-1, j:j+patch_size-1) = ...
weight(i:i+patch_size-1, j:j+patch_size-1) + 1;
end
end
% 平均重叠区域
denoised_img = denoised_img ./ weight;
end
6. 性能优化技巧
内存优化版本
function [L, S] = rpca_memory_efficient(M, lambda, max_iter)
% 内存高效的RPCA实现
[m, n] = size(M);
% 使用随机SVD加速计算
opts.tol = 1e-4;
opts.maxit = 50;
L = zeros(m, n);
S = zeros(m, n);
Y = zeros(m, n);
mu = 1.25 / norm(M, 2);
for iter = 1:max_iter
% 使用随机SVD计算低秩近似
[U, s, V] = rsvd(M - S + Y/mu, min(m,n), opts);
s = max(s - 1/mu, 0);
L = U * diag(s) * V';
% 更新稀疏成分
S = soft_threshold(M - L + Y/mu, lambda/mu);
% 更新乘子
Y = Y + mu * (M - L - S);
% 简单的收敛检查
if norm(M - L - S, 'fro') < 1e-6 * norm(M, 'fro')
break;
end
end
end
function X = soft_threshold(X, tau)
% 软阈值函数
X = sign(X) .* max(abs(X) - tau, 0);
end
function [U, S, V] = rsvd(A, k, opts)
% 随机SVD计算
if nargin < 3
opts.tol = 1e-4;
opts.maxit = 50;
end
[m, n] = size(A);
Omega = randn(n, k);
Y = A * Omega;
[Q, ~] = qr(Y, 0);
B = Q' * A;
[U_tilde, S, V] = svd(B, 'econ');
U = Q * U_tilde;
U = U(:, 1:k);
S = diag(S);
S = S(1:k);
V = V(:, 1:k);
end
关键参数选择建议
- λ参数：通常设为
1/sqrt(max(m,n)) - μ参数：通常设为
1.25/norm(M,2) - 收敛容差：
1e-6到1e-8 - 最大迭代次数：100-1000，取决于图像大小
这种方法特别适合处理包含稀疏大噪声的图像，能够有效分离图像的结构信息和噪声成分。
