✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,代码获取、论文复现及科研仿真合作可私信。
🍎个人主页:Matlab科研工作室
🍊个人信条:格物致知。
更多Matlab完整代码及仿真定制内容点击👇
🔥 内容介绍
图像去噪是图像处理领域中一项重要的基础任务,旨在从图像中去除噪声,提高图像质量。近年来,随着深度学习技术的快速发展,基于深度学习的图像去噪方法取得了显著的成果。然而,深度学习方法通常需要大量的训练数据,并且对模型结构和参数设置较为敏感。
半盲图像去噪是指只知道部分图像信息的情况下进行图像去噪,例如只知道图像的梯度信息或边缘信息。半盲图像去噪在许多实际应用中都非常重要,例如图像修复、图像增强和图像压缩等。
基于TV算法的半盲图像去噪
全变分 (Total Variation, TV) 算法是一种经典的图像去噪方法,它利用图像梯度的稀疏性来恢复图像。TV算法的数学模型如下:
min_u ||∇u||_1 + λ||u - f||_2^2
其中,u 是去噪后的图像,f 是带噪声的图像,∇u 是 u 的梯度,λ 是正则化参数。
对于半盲图像去噪,我们只知道图像的梯度信息或边缘信息,因此需要对 TV 算法进行修改。一种常用的方法是将梯度信息或边缘信息作为约束条件添加到 TV 算法中。例如,如果我们只知道图像的梯度信息,则可以将 TV 算法修改为:
min_u ||∇u||_1 + λ||∇u - g||_2^2
其中,g 是图像的梯度信息.
基于 TV 算法的半盲图像去噪方法是一种有效且实用的图像去噪方法,它能够在只知道部分图像信息的情况下有效地去除噪声,提高图像质量。该方法在图像修复、图像增强和图像压缩等实际应用中具有重要的应用价值。
📣 部分代码
%% Experiment - Deblurring with Total-Variation prior: Laplace point spread function as kernel.
% -------------------------------------------------------------------------%
% -------------------------------------------------------------------------%
clear all;clc;
fix_sigma = 0;
op.fix_b = fix_b;
op.fix_sigma = fix_sigma;
%% RUN EXPERIMENTS
for b = [0.3]
op.b = b;
if fix_b
fprintf('\n\n\t\t Fix b \n\n');
op.b_init = b;
end
for snr = [30]
op.BSNR = snr;
fprintf('\n\n===========<< SNR = %d \t|\t B = %d>>============\n\n', snr, b);
for i_im = [8]
filename=testImg{i_im};
[~,fix,~] = fileparts(char(filename));
name = fix;
%% Generate a noisy and blurred observation vector 'y'
% original image 'x'
fprintf(['File:' filename '\n']);
x = double(imread(strcat('images/', filename)));
op.x = x;
dimX = numel(x);
im_size = size(x);
% kernel function
H_FFT = @(b) laplace_psf(im_size, psf_size, b);
% conjugate of the PSF
HC_FFT = @(b) conj(H_FFT(b));
% differential w.r.t b
diff_fftkernel_b = @(b) diff_laplace_b(im_size, psf_size, b);
%% psf and its gradients
A = @(x, b) real(ifft2(H_FFT(b) .* fft2(x)));
AT = @(x, b) real(ifft2(HC_FFT(b) .* fft2(x)));
diff_A_b = @(x, b) real(ifft2(diff_fftkernel_b(b) .* fft2(x)));
%% Max eigen value of A through iterative method
evMax=max_eigenval_Laplace(A,AT, 1, im_size, 1e-4, 1e4, 0);
fprintf('evmax=%d\n',evMax);
%% observation y
Ax = real(A(x, b));
sigma = norm(Ax-mean(mean(Ax)),'fro')/sqrt(dimX*10^(op.BSNR/10));
fprintf('Sigma^2 = %d\t\t snr:%d\n\n', sigma^2, op.BSNR);
op.sigma = sigma;
sigma_min= norm(Ax-mean(mean(Ax)),'fro')/sqrt(dimX*10^(op.BSNR_min/10));
sigma_max = norm(Ax-mean(mean(Ax)),'fro')/sqrt(dimX*10^(op.BSNR_max/10));
op.sigma_init = (sigma_min^2 + sigma_max^2) / 2;
if fix_sigma
op.sigma_init = sigma^2;
end
op.sigma_min = sigma_min^2;
op.sigma_max = sigma_max^2;
y = Ax + sigma * randn(size(Ax));
op.X0 = y;
%% Experiment setup: functions related to Bayesian model
% Likelihood and gradients
op.f = @(x, b, sigma2) (norm(y-A(x, b),'fro')^2)/(2*sigma2);
op.gradF = @(x, b, sigma2) real(AT(A(x,b) - y, b)/sigma2);
op.grad_b = @(x, b, sigma2) real(sum(sum(diff_A_b(x,b) .* (A(x, b) - y))) / sigma2);
op.gradF_sigma = @(x, b, sigma2) ((norm(y - A(x, b) ,'fro')^2)/(2*sigma2^2) - dimX/(2*sigma2)) ;
Lf = @ (sigma2) evMax^2 / sigma2;
op.Lf = max(Lf(sigma_min^2), Lf(sigma_max^2));
% Set algorithm parameters that depend on Lf
op.lambda= min(5/op.Lf,op.lambdaMax);%smoothing parameter for MYULA sampler
op.gamma_max = 1 / (op.Lf+(1/op.lambda));%Lf is the lipschitz constant of the gradient of the smooth part
op.gamma= 10 * op.gammaFrac * op.gamma_max;%discretisation step MYULA
% Regulariser % TV norm
op.g = @(x) TVnorm(x); %g(x) for TV regularizer
chambolleit = 25;
% Proximal operator of g(x)
op.proxG = @(x,lambda,theta) chambolle_prox_TV_stop(x,'lambda',lambda*theta,'maxiter',chambolleit);
Psi = @(x,th) op.proxG(x, 1, th); % define this format for SALSA solver
% We use this scalar summary to monitor convergence
op.logPi = @(x,theta, b, sigma2) -op.f(x, b, sigma2) -theta*op.g(x);
%% Run SAPG Algorithm 1 to compute theta_EB
tic;
[theta_EB, b_EB, sigma2_EB, results] = SAPG_algorithm_laplace(y, op);
SAPG_time = toc;
fprintf('\n\n\t\t=========>>> SAPG Elapsed: %d <<<=========\n\n', SAPG_time);
results.SAPG_time = SAPG_time;
results.theta_EB = theta_EB;
results.b_EB = b_EB;
results.sigma2_EB = sigma2_EB;
results.b = b;
sigma_EB = sigma2_EB;
results.sigma = sigma;
last_b = results.last_b;
results.op = op;
results.x=x;
results.y=y;
fprintf('Theta EB: %d\t|\t B EB: %d\t|\t Last b: %d\t|\t Sigma EB: %d', theta_EB, b_EB, last_b, sigma_EB);
%% Solve MYULA problem with theta_EB, b_EB, mu_EB and sigma_EB
% b_EB, sigma_EB, theta_EB,
tic;
global calls;
calls = 0;
outeriters = 500;
A1 = @(x) A(x, b_EB);
AT1 = @(x) AT(x, b_EB);
A1 = @(x) callcounter(A1,x);
AT1 = @(x) callcounter(AT1,x);
inneriters = 1;
tol = 1e-5; %taken from salsa demo.
mu_ = theta_EB/10;
muinv = 1/mu_;
filter_FFT = 1./(abs(H_FFT(b_EB)).^2 + mu_);
invLS = @(x) real(ifft2(filter_FFT.*fft2( x)));
invLS = @(xw) callcounter(invLS,xw);
fprintf('Running SALSA solver...\n')
% SALSA
[xMAP, ~, ~, ~, ~, ~, ~] = ...
SALSA_v2(y, A1, theta_EB*sigma_EB,...
'MU', mu_, ...
'AT', AT1, ...
'StopCriterion', 1, ...
'True_x', x, ...
'ToleranceA', tol, ...
'MAXITERA', outeriters, ...
'Psi', Psi, ...
'Phi', op.g, ...
'TVINITIALIZATION', 1, ...
'TViters', 10, ...
'LS', invLS, ...
'VERBOSE', 0);
%Compute MSE with xMAP
mse_b = 10*log10(norm(x-xMAP,'fro')^2 /dimX);
results.mse_b=mse_b;
results.xMAP=xMAP;
results.snr_true_reconst_image_with_estimate_para = snr_func(x, xMAP);
[ssimval, ~] = ssim(x, xMAP);
results.ssim_true_reconst_image_with_estimate_para = ssimval;
% Sigma^2
figSigma=figure;
sig_true = results.sigma^2 * ones(results.last_samp,1);
plot(results.sigmas(1:results.last_samp),'b','LineWidth',1.5,'MarkerSize',8);hold on;
plot(sig_true, 'r--','LineWidth',1.5,'MarkerSize',8);
xlim([1 inf]); grid on; xlabel('Iteration (n)');ylabel('$\sigma^2_n$', 'Interpreter', 'latex');title(strcat('$\sigma^2_{EB} = $',num2str(sigma_EB)), 'Interpreter', 'latex');
legend('$\sigma^2_n$', '$\sigma^2_{true}$', 'Interpreter', 'latex');hold off;
% Theta
figTheta=figure;
plot(results.thetas(1:results.last_samp),'b','LineWidth',1.5,'MarkerSize',8);
xlim([1 inf]); grid on; xlabel('Iteration (n)');ylabel('\theta_n');title(strcat('$\theta_{EB} = $',num2str(theta_EB)), 'Interpreter', 'latex');
% b
true_b = results.b * ones(results.last_samp,1);
figB=figure;
plot(results.bs(1:results.last_samp),'b','LineWidth',1.5,'MarkerSize',8);hold on;
plot(true_b, 'r--','LineWidth',1.5,'MarkerSize',8);hold off;
legend('$b_n$', '$b_{true}$', 'Interpreter', 'latex');
xlim([1 inf]); grid on; xlabel('Iteration (n)');ylabel('b_n');title(strcat('$b_{EB}$ = ',num2str(b_EB)), 'Interpreter', 'latex');
% images
figtrue = figure;
imagesc(x), colormap gray, axis off, axis equal;
figmean.Units = 'inches';
figmean.OuterPosition = [0.25 0.25 10 10];
title("x")
figtrue = figure;
imagesc(y), colormap gray, axis off, axis equal;
figmean.Units = 'inches';
figmean.OuterPosition = [0.25 0.25 10 10];
title("y")
figMap = figure;
imagesc(xMAP), colormap gray, axis off, axis equal;
figMap.Units = 'inches';
figMap.OuterPosition = [0.25 0.25 10 10];
title("x MAP")
end
end
end
⛳️ 运行结果
🔗 参考文献
[1] 王满.针对全变分图像去噪的半光滑牛顿法研究[D].昆明理工大学,2017.DOI:CNKI:CDMD:2.1017.222587.
[2] 来沛剑.基于全变分的含噪图像盲复原算法研究[M].[2024-04-29].
[3] 陈晓钢.基于统计分析的图像去模糊与图像去噪新方法研究[D].上海交通大学[2024-04-29].
🎈 部分理论引用网络文献,若有侵权联系博主删除
🎁 关注我领取海量matlab电子书和数学建模资料
👇 私信完整代码和数据获取及论文数模仿真定制
1 各类智能优化算法改进及应用
生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化、背包问题、 风电场布局、时隙分配优化、 最佳分布式发电单元分配、多阶段管道维修、 工厂-中心-需求点三级选址问题、 应急生活物质配送中心选址、 基站选址、 道路灯柱布置、 枢纽节点部署、 输电线路台风监测装置、 集装箱船配载优化、 机组优化、 投资优化组合、云服务器组合优化、 天线线性阵列分布优化、CVRP问题、VRPPD问题、多中心VRP问题、多层网络的VRP问题、多中心多车型的VRP问题、 动态VRP问题、双层车辆路径规划(2E-VRP)、充电车辆路径规划(EVRP)、油电混合车辆路径规划、混合流水车间问题、 订单拆分调度问题、 公交车的调度排班优化问题、航班摆渡车辆调度问题、选址路径规划问题
2 机器学习和深度学习方面
2.1 bp时序、回归预测和分类
2.2 ENS声神经网络时序、回归预测和分类
2.3 SVM/CNN-SVM/LSSVM/RVM支持向量机系列时序、回归预测和分类
2.4 CNN/TCN卷积神经网络系列时序、回归预测和分类
2.5 ELM/KELM/RELM/DELM极限学习机系列时序、回归预测和分类
2.6 GRU/Bi-GRU/CNN-GRU/CNN-BiGRU门控神经网络时序、回归预测和分类
2.7 ELMAN递归神经网络时序、回归\预测和分类
2.8 LSTM/BiLSTM/CNN-LSTM/CNN-BiLSTM/长短记忆神经网络系列时序、回归预测和分类
2.9 RBF径向基神经网络时序、回归预测和分类