【图像去噪】基于TV算法实现半盲图像去噪附matlab代码

 ✅作者简介：热爱科研的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 EXPERIMENTSfor 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  endend​​

⛳️ 运行结果

🔗 参考文献

[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径向基神经网络时序、回归预测和分类

2.10 DBN深度置信网络时序、回归预测和分类

2.11 FNN模糊神经网络时序、回归预测

2.12 RF随机森林时序、回归预测和分类

2.13 BLS宽度学习时序、回归预测和分类

2.14 PNN脉冲神经网络分类

2.15 模糊小波神经网络预测和分类

2.16 时序、回归预测和分类

2.17 时序、回归预测预测和分类

2.18 XGBOOST集成学习时序、回归预测预测和分类

方向涵盖风电预测、光伏预测、电池寿命预测、辐射源识别、交通流预测、负荷预测、股价预测、PM2.5浓度预测、电池健康状态预测、用电量预测、水体光学参数反演、NLOS信号识别、地铁停车精准预测、变压器故障诊断

2.图像处理方面

图像识别、图像分割、图像检测、图像隐藏、图像配准、图像拼接、图像融合、图像增强、图像压缩感知

3 路径规划方面

旅行商问题（TSP）、车辆路径问题（VRP、MVRP、CVRP、VRPTW等）、无人机三维路径规划、无人机协同、无人机编队、机器人路径规划、栅格地图路径规划、多式联运运输问题、 充电车辆路径规划（EVRP）、 双层车辆路径规划（2E-VRP）、 油电混合车辆路径规划、 船舶航迹规划、 全路径规划规划、 仓储巡逻

4 无人机应用方面

无人机路径规划、无人机控制、无人机编队、无人机协同、无人机任务分配、无人机安全通信轨迹在线优化、车辆协同无人机路径规划

5 无线传感器定位及布局方面

传感器部署优化、通信协议优化、路由优化、目标定位优化、Dv-Hop定位优化、Leach协议优化、WSN覆盖优化、组播优化、RSSI定位优化

6 信号处理方面

信号识别、信号加密、信号去噪、信号增强、雷达信号处理、信号水印嵌入提取、肌电信号、脑电信号、信号配时优化

7 电力系统方面

微电网优化、无功优化、配电网重构、储能配置、有序充电

8 元胞自动机方面

交通流 人群疏散 病毒扩散 晶体生长 金属腐蚀

9 雷达方面

卡尔曼滤波跟踪、航迹关联、航迹融合