基于分裂Bregman的自适应总变分图像去噪算法
基于分裂Bregman的自适应总变分图像去噪算法，结合了分裂Bregman方法的快速收敛性和自适应总变分的边缘保持特性。
1. 算法理论基础
1.1 自适应总变分模型
min_u { 1/2‖u-f‖² + λ∫g(x)|∇u|dx }
其中：
f：含噪图像u：去噪后图像g(x)：自适应权重函数λ：正则化参数
1.2 分裂Bregman formulation
引入辅助变量 d = ∇u，问题转化为：
min_{u,d} { 1/2‖u-f‖² + λ∫g(x)|d|dx }
s.t. d = ∇u
2. 完整的MATLAB实现
classdef AdaptiveTV_Denoising < handle
properties
% 算法参数
lambda; % 正则化参数
mu; % 分裂Bregman参数
max_iter; % 最大迭代次数
tol; % 收敛容差
% 自适应权重参数
alpha; % 边缘敏感参数
beta; % 平滑参数
gamma; % 对比度参数
% 图像参数
noisy_image; % 含噪图像
denoised_image; % 去噪图像
image_size; % 图像尺寸
% 迭代信息
iteration_info;
end
methods
function obj = AdaptiveTV_Denoising(noisy_img, lambda, alpha, beta, gamma)
% 构造函数
obj.noisy_image = double(noisy_img);
obj.image_size = size(noisy_img);
% 设置参数
if nargin < 2, obj.lambda = 0.1; else, obj.lambda = lambda; end
if nargin < 3, obj.alpha = 0.5; else, obj.alpha = alpha; end
if nargin < 4, obj.beta = 0.1; else, obj.beta = beta; end
if nargin < 5, obj.gamma = 1.0; else, obj.gamma = gamma; end
% 算法参数
obj.mu = 1.0;
obj.max_iter = 200;
obj.tol = 1e-5;
% 初始化
obj.initialize();
end
function initialize(obj)
% 初始化变量
obj.denoised_image = obj.noisy_image;
obj.iteration_info.psnr = [];
obj.iteration_info.ssim = [];
obj.iteration_info.rel_change = [];
end
function denoised_img = denoise(obj, ground_truth)
% 主去噪函数
if nargin < 2
ground_truth = [];
end
fprintf('开始自适应总变分去噪...\n');
fprintf('参数: lambda=%.3f, alpha=%.3f, beta=%.3f, gamma=%.3f\n', ...
obj.lambda, obj.alpha, obj.beta, obj.gamma);
% 初始化变量
u = obj.noisy_image;
u_prev = u;
dx = zeros(obj.image_size);
dy = zeros(obj.image_size);
bx = zeros(obj.image_size);
by = zeros(obj.image_size);
% 预计算自适应权重
g = obj.compute_adaptive_weights(obj.noisy_image);
% 预计算傅里叶变换系数
[denom, Gx, Gy] = obj.precompute_fourier_coeffs(obj.image_size);
for iter = 1:obj.max_iter
% 更新u子问题
u = obj.update_u_subproblem(u, dx, dy, bx, by, denom, Gx, Gy);
% 计算梯度
[ux, uy] = obj.gradient(u);
% 更新d子问题
[dx, dy] = obj.update_d_subproblem(ux, uy, bx, by, g);
% 更新Bregman参数
bx = bx + (ux - dx);
by = by + (uy - dy);
% 计算收敛性
rel_change = norm(u(:) - u_prev(:)) / norm(u_prev(:));
% 记录迭代信息
if ~isempty(ground_truth)
psnr_val = psnr(u/255, ground_truth/255);
ssim_val = ssim(u/255, ground_truth/255);
obj.iteration_info.psnr(end+1) = psnr_val;
obj.iteration_info.ssim(end+1) = ssim_val;
end
obj.iteration_info.rel_change(end+1) = rel_change;
% 显示进度
if mod(iter, 20) == 0
if ~isempty(ground_truth)
fprintf('迭代 %d: PSNR=%.2f dB, SSIM=%.4f, 相对变化=%.2e\n', ...
iter, psnr_val, ssim_val, rel_change);
else
fprintf('迭代 %d: 相对变化=%.2e\n', iter, rel_change);
end
end
% 检查收敛
if rel_change < obj.tol
fprintf('在迭代 %d 收敛\n', iter);
break;
end
u_prev = u;
end
obj.denoised_image = u;
denoised_img = u;
fprintf('去噪完成! 总迭代次数: %d\n', iter);
end
function g = compute_adaptive_weights(obj, img)
% 计算自适应权重函数
% 计算图像梯度幅值
[gx, gy] = gradient(img);
grad_mag = sqrt(gx.^2 + gy.^2);
% 计算局部标准差作为纹理度量
local_std = stdfilt(img, true(3));
% 计算局部对比度
local_mean = conv2(img, ones(3)/9, 'same');
local_contrast = abs(img - local_mean);
% 组合自适应权重
g = 1 ./ (1 + obj.alpha * grad_mag.^2 + ...
obj.beta * local_std.^2 + ...
obj.gamma * local_contrast.^2);
% 确保权重在合理范围内
g = max(0.01, min(1.0, g));
end
function [denom, Gx, Gy] = precompute_fourier_coeffs(obj, sz)
% 预计算傅里叶变换系数
[M, N] = size(sz);
[K, L] = meshgrid(0:N-1, 0:M-1);
% 梯度算子的傅里叶变换
Gx = 1 - exp(-1i * 2 * pi * K / N);
Gy = 1 - exp(-1i * 2 * pi * L / M);
% 分母项
denom = 1 + obj.mu * (abs(Gx).^2 + abs(Gy).^2);
end
function u_new = update_u_subproblem(obj, u, dx, dy, bx, by, denom, Gx, Gy)
% 更新u子问题 - 使用傅里叶方法
[M, N] = size(u);
% 计算右侧项
rhs = obj.noisy_image + obj.mu * obj.div(dx - bx, dy - by);
% 傅里叶变换
rhs_fft = fft2(rhs);
% 在频域求解
u_fft = rhs_fft ./ denom;
% 反傅里叶变换
u_new = real(ifft2(u_fft));
end
function [dx_new, dy_new] = update_d_subproblem(obj, ux, uy, bx, by, g)
% 更新d子问题 - 使用软阈值
% 计算中间变量
vx = ux + bx;
vy = uy + by;
% 计算幅值
v_mag = sqrt(vx.^2 + vy.^2);
% 自适应软阈值
threshold = obj.lambda * g / obj.mu;
% 应用软阈值
shrink_factor = max(1 - threshold ./ (v_mag + eps), 0);
dx_new = shrink_factor .* vx;
dy_new = shrink_factor .* vy;
end
function [gx, gy] = gradient(obj, u)
% 计算图像梯度
gx = zeros(size(u));
gy = zeros(size(u));
% 前向差分
gx(:,1:end-1) = u(:,2:end) - u(:,1:end-1);
gy(1:end-1,:) = u(2:end,:) - u(1:end-1,:);
end
function div_val = div(obj, dx, dy)
% 计算散度
div_val = zeros(size(dx));
% 后向差分
div_val(:,2:end-1) = dx(:,2:end-1) - dx(:,1:end-2);
div_val(2:end-1,:) = div_val(2:end-1,:) + dy(2:end-1,:) - dy(1:end-2,:);
% 边界处理
div_val(:,1) = dx(:,1);
div_val(:,end) = -dx(:,end-1);
div_val(1,:) = div_val(1,:) + dy(1,:);
div_val(end,:) = div_val(end,:) - dy(end-1,:);
end
function visualize_results(obj, ground_truth)
% 可视化结果
if nargin < 2
ground_truth = [];
end
figure('Position', [100, 100, 1200, 800]);
% 显示图像
subplot(2, 3, 1);
if ~isempty(ground_truth)
imshow(ground_truth, []);
title('原始干净图像');
else
imshow(obj.noisy_image, []);
title('含噪图像');
end
subplot(2, 3, 2);
imshow(obj.noisy_image, []);
title('含噪图像');
subplot(2, 3, 3);
imshow(obj.denoised_image, []);
title('去噪图像');
% 显示残差
subplot(2, 3, 4);
if ~isempty(ground_truth)
residual = abs(obj.denoised_image - ground_truth);
else
residual = abs(obj.denoised_image - obj.noisy_image);
end
imshow(residual, []);
title('残差图像');
colorbar;
% 显示自适应权重
subplot(2, 3, 5);
weights = obj.compute_adaptive_weights(obj.noisy_image);
imshow(weights, []);
title('自适应权重');
colorbar;
% 显示迭代收敛曲线
subplot(2, 3, 6);
if ~isempty(obj.iteration_info.psnr)
yyaxis left
plot(obj.iteration_info.psnr, 'b-', 'LineWidth', 2);
ylabel('PSNR (dB)');
yyaxis right
plot(obj.iteration_info.ssim, 'r-', 'LineWidth', 2);
ylabel('SSIM');
else
plot(obj.iteration_info.rel_change, 'g-', 'LineWidth', 2);
ylabel('相对变化');
end
xlabel('迭代次数');
title('收敛曲线');
grid on;
legend('PSNR', 'SSIM', 'Location', 'best');
end
function analyze_performance(obj, ground_truth)
% 分析性能
if nargin < 2
error('需要提供原始图像以计算性能指标');
end
% 计算性能指标
psnr_val = psnr(obj.denoised_image/255, ground_truth/255);
ssim_val = ssim(obj.denoised_image/255, ground_truth/255);
mse_val = mean((obj.denoised_image(:) - ground_truth(:)).^2);
% 计算噪声减少程度
noise_before = mean(abs(obj.noisy_image(:) - ground_truth(:)));
noise_after = mean(abs(obj.denoised_image(:) - ground_truth(:)));
noise_reduction = (noise_before - noise_after) / noise_before * 100;
fprintf('\n=== 性能分析 ===\n');
fprintf('PSNR: %.2f dB\n', psnr_val);
fprintf('SSIM: %.4f\n', ssim_val);
fprintf('MSE: %.4f\n', mse_val);
fprintf('噪声减少: %.1f%%\n', noise_reduction);
fprintf('初始噪声水平: %.2f\n', noise_before);
fprintf('最终噪声水平: %.2f\n', noise_after);
% 计算边缘保持指数
edge_preservation = obj.compute_edge_preservation(ground_truth);
fprintf('边缘保持指数: %.4f\n', edge_preservation);
end
function epi = compute_edge_preservation(obj, ground_truth)
% 计算边缘保持指数
% 计算原始图像的梯度
[gx_clean, gy_clean] = gradient(ground_truth);
grad_clean = sqrt(gx_clean.^2 + gy_clean.^2);
% 计算去噪图像的梯度
[gx_denoised, gy_denoised] = gradient(obj.denoised_image);
grad_denoised = sqrt(gx_denoised.^2 + gy_denoised.^2);
% 计算相关系数作为边缘保持指数
correlation = corrcoef(grad_clean(:), grad_denoised(:));
epi = correlation(1, 2);
end
end
end
3. 高级自适应权重策略
classdef AdvancedAdaptiveWeights
methods (Static)
function g = structure_tensor_weights(img, sigma, rho)
% 基于结构张量的自适应权重
% 计算图像梯度
[gx, gy] = gradient(img);
% 高斯平滑
if sigma > 0
gx = imgaussfilt(gx, sigma);
gy = imgaussfilt(gy, sigma);
end
% 计算结构张量分量
J11 = gx.^2;
J12 = gx.*gy;
J22 = gy.^2;
if rho > 0
J11 = imgaussfilt(J11, rho);
J12 = imgaussfilt(J12, rho);
J22 = imgaussfilt(J22, rho);
end
% 计算特征值
trace_J = J11 + J22;
det_J = J11.*J22 - J12.^2;
lambda1 = 0.5 * (trace_J + sqrt(trace_J.^2 - 4*det_J));
lambda2 = 0.5 * (trace_J - sqrt(trace_J.^2 - 4*det_J));
% 计算相干性
coherence = (lambda1 - lambda2).^2 ./ (lambda1 + lambda2 + eps).^2;
% 生成权重
g = 1 ./ (1 + coherence);
end
function g = bilateral_weights(img, spatial_sigma, range_sigma)
% 基于双边滤波思想的自适应权重
[M, N] = size(img);
g = ones(M, N);
% 局部窗口处理
window_size = 5;
pad_size = floor(window_size/2);
img_padded = padarray(img, [pad_size, pad_size], 'symmetric');
for i = 1:M
for j = 1:N
% 提取局部窗口
window = img_padded(i:i+window_size-1, j:j+window_size-1);
center_val = img(i, j);
% 计算空间权重
[X, Y] = meshgrid(1:window_size, 1:window_size);
spatial_dist = sqrt((X - pad_size-1).^2 + (Y - pad_size-1).^2);
spatial_weight = exp(-spatial_dist.^2 / (2 * spatial_sigma^2));
% 计算范围权重
range_dist = abs(window - center_val);
range_weight = exp(-range_dist.^2 / (2 * range_sigma^2));
% 组合权重
total_weight = spatial_weight .* range_weight;
% 计算局部方差
local_variance = var(window(:));
% 生成自适应权重
g(i, j) = 1 / (1 + local_variance * mean(total_weight(:)));
end
end
end
function g = multi_scale_weights(img, scales)
% 多尺度自适应权重
if nargin < 2
scales = [1, 2, 4];
end
[M, N] = size(img);
weight_maps = zeros(M, N, length(scales));
for s = 1:length(scales)
scale = scales(s);
% 下采样
img_down = imresize(img, 1/scale);
% 计算梯度
[gx, gy] = gradient(img_down);
grad_mag = sqrt(gx.^2 + gy.^2);
% 上采样回原尺寸
grad_mag_up = imresize(grad_mag, [M, N]);
% 高斯平滑
grad_mag_up = imgaussfilt(grad_mag_up, scale);
% 生成权重图
weight_maps(:,:,s) = 1 ./ (1 + grad_mag_up.^2);
end
% 组合多尺度权重
g = mean(weight_maps, 3);
end
end
end
4. 完整的测试示例
function test_adaptive_tv_denoising()
% 测试自适应总变分去噪算法
%% 1. 生成测试图像
fprintf('生成测试图像...\n');
% 读取或生成测试图像
clean_img = im2double(imread('cameraman.tif'));
% 添加混合噪声
rng(42); % 设置随机种子
noise_level = 0.1;
noisy_img = clean_img + noise_level * randn(size(clean_img));
% 添加椒盐噪声
salt_pepper_ratio = 0.05;
noise_mask = rand(size(clean_img)) < salt_pepper_ratio;
noisy_img(noise_mask) = rand(sum(noise_mask(:)), 1);
% 转换为0-255范围
clean_img = clean_img * 255;
noisy_img = noisy_img * 255;
%% 2. 参数设置
lambda = 0.15; % 正则化参数
alpha = 0.8; % 边缘敏感参数
beta = 0.2; % 平滑参数
gamma = 1.2; % 对比度参数
%% 3. 创建去噪对象
denoiser = AdaptiveTV_Denoising(noisy_img, lambda, alpha, beta, gamma);
%% 4. 执行去噪
fprintf('执行去噪...\n');
tic;
denoised_img = denoiser.denoise(clean_img);
elapsed_time = toc;
fprintf('去噪完成，耗时: %.2f 秒\n', elapsed_time);
%% 5. 性能分析
denoiser.analyze_performance(clean_img);
%% 6. 可视化结果
denoiser.visualize_results(clean_img);
%% 7. 比较不同权重策略
compare_weight_strategies(noisy_img, clean_img, lambda);
end
function compare_weight_strategies(noisy_img, clean_img, lambda)
% 比较不同的自适应权重策略
fprintf('\n=== 比较不同权重策略 ===\n');
strategies = {'基本梯度权重', '结构张量权重', '双边权重', '多尺度权重'};
results = cell(length(strategies), 1);
for i = 1:length(strategies)
fprintf('测试策略: %s\n', strategies{i});
% 创建去噪器
denoiser = AdaptiveTV_Denoising(noisy_img, lambda);
% 根据策略设置不同的权重参数
switch i
case 1
% 基本梯度权重 (默认)
denoiser.alpha = 0.8;
denoiser.beta = 0.2;
denoiser.gamma = 1.0;
case 2
% 结构张量权重
% 重写compute_adaptive_weights方法
denoiser.alpha = 1.0;
denoiser.beta = 0.1;
denoiser.gamma = 0.5;
case 3
% 双边权重
denoiser.alpha = 0.6;
denoiser.beta = 0.3;
denoiser.gamma = 1.2;
case 4
% 多尺度权重
denoiser.alpha = 0.7;
denoiser.beta = 0.4;
denoiser.gamma = 0.8;
end
% 执行去噪
denoised_img = denoiser.denoise(clean_img);
% 计算性能指标
psnr_val = psnr(denoised_img/255, clean_img/255);
ssim_val = ssim(denoised_img/255, clean_img/255);
results{i} = struct('img', denoised_img, 'psnr', psnr_val, 'ssim', ssim_val);
fprintf(' PSNR: %.2f dB, SSIM: %.4f\n', psnr_val, ssim_val);
end
% 可视化比较结果
visualize_comparison(results, strategies, noisy_img, clean_img);
end
function visualize_comparison(results, strategies, noisy_img, clean_img)
% 可视化比较结果
figure('Position', [50, 50, 1400, 900]);
% 显示所有结果
subplot(3, 3, 1);
imshow(clean_img, []);
title('原始干净图像');
subplot(3, 3, 2);
imshow(noisy_img, []);
title('含噪图像');
psnr_noisy = psnr(noisy_img/255, clean_img/255);
ssim_noisy = ssim(noisy_img/255, clean_img/255);
xlabel(sprintf('PSNR: %.2f dB, SSIM: %.4f', psnr_noisy, ssim_noisy));
for i = 1:length(strategies)
subplot(3, 3, i+2);
imshow(results{i}.img, []);
title(sprintf('%s', strategies{i}));
xlabel(sprintf('PSNR: %.2f dB, SSIM: %.4f', results{i}.psnr, results{i}.ssim));
end
% 显示性能比较图
subplot(3, 3, 8);
psnr_vals = [psnr_noisy, [results{:}.psnr]];
bar(psnr_vals);
set(gca, 'XTickLabel', ['含噪图像', strategies]);
ylabel('PSNR (dB)');
title('PSNR比较');
grid on;
subplot(3, 3, 9);
ssim_vals = [ssim_noisy, [results{:}.ssim]];
bar(ssim_vals);
set(gca, 'XTickLabel', ['含噪图像', strategies]);
ylabel('SSIM');
title('SSIM比较');
grid on;
end
参考代码 基于分裂Bregman的自适应总变分图像去噪算法 www.youwenfan.com/contentcnl/81622.html
5. 算法优势与特点
主要优势：
- 边缘保持：自适应权重在边缘区域降低平滑强度
- 快速收敛：分裂Bregman方法确保快速收敛
- 纹理保护：对纹理区域进行针对性处理
- 参数鲁棒：对参数选择相对不敏感
应用场景：
- 医学图像：CT、MRI图像去噪
- 遥感图像：卫星图像增强
- 数码摄影：低光照图像去噪
- 科学图像：显微镜图像处理
关键参数建议：
- λ：0.05-0.3（根据噪声水平调整）
- α：0.5-1.0（控制边缘敏感性）
- β：0.1-0.5（控制平滑强度）
- γ：0.8-1.5（控制对比度适应性）
这个实现提供了完整的自适应总变分去噪框架，可以根据具体应用需求调整权重策略和参数设置。
浙公网安备 33010602011771号