【图像去噪】基于柯西近端分裂 (CPS) 算法实现图像去噪附MATLAB源代码

1 简介

In this paper, we propose a proximal splitting methodology with a non-convex penalty function based on the heavy-tailed Cauchy distribution. We first suggest a closed-form expression for calculating the proximal operator of the Cauchy prior, which then makes it applicable in generic proximal splitting algorithms. We further derive the condition required for guaranteed convergence to the global minimum in optimisation problems involving the Cauchy based penalty function. Setting the system parameters by satisfying the proposed condition ensures convergence even though the overall cost function is non-convex, when minimisation is performed via a proximal splitting algorithm. The proposed method based on Cauchy regularisation is evaluated by solving generic signal processing examples, i.e. 1D signal denoising in the frequency domain, two image reconstruction tasks including de-blurring and denoising, and error recovery in a multiple-antenna communication system. We experimentally verify the proposed convergence conditions for various cases, and show the effectiveness of the proposed Cauchy based non-convex penalty function over state-of-the-art penalty functions such as $L_1$ and total variation ( $TV$ ) norms.

2 部分代码

%% Deblurring via Cauchy proximal splitting algorithm% Y = HX + N% Y is the blurred and noisy image% X is the clear image (object of interest)% H is the point spread function (blurring operation)% N is the additive zero-meam Gaussian noise with variance \sigma^2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Some Important Variables%       ** X: Noise-free Image.%%       ** filterSize: The size of the blurring filter.%%       ** BSNRdb: Blurred-SNR value in decibels.%%       ** maxIter: Maximum number of FB iterations.%%       ** Y: Noisy and blurred image.%%       ** mu: CPS step size%%       ** gamma: Cauchy scale parameter%%       ** errorCriterion: The stopping criterion for the FB algorithm.%%       ** x_hat: The reconstructed image.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LICENSE% % This program is free software: you can redistribute it and/or modify% it under the terms of the GNU General Public License as published by% the Free Software Foundation, either version 3 of the License, or% (at your option) any later version.% % This program is distributed in the hope that it will be useful,% but WITHOUT ANY WARRANTY; without even the implied warranty of% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the% GNU General Public License for more details.% % You should have received a copy of the GNU General Public License% along with this program.  If not, see <https://www.gnu.org/licenses/>.% % Copyright (C) Oktay Karakus,PhD % University of Bristol, UK% o.karakus@bristol.ac.uk% April 2020%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% REFERENCE%% [1] O Karakus, P Mayo, and A Achim. "Convergence Guarantees for %     Non-Convex Optimisation with Cauchy-Based Penalties"%       arXiv preprint.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clearvarsclose allclcimName = 'cameraman';X = double(imread([imName '.tif']));filterSize = 5;h = fspecial('gaussian',[filterSize filterSize], 1);blurred = imfilter(X, h, 'circular');BSNRdb = 40; % Blurred-signal to noise ratio in decibelssigma = norm(blurred-mean(mean(blurred)),'fro')/sqrt(size(blurred, 1)*size(blurred, 2)*10^(BSNRdb/10)); % noise std calculationY = blurred + sigma*randn(size(blurred)); % The noisy and blurred image Y.K  = psf2otf(h, size(X));KC = conj(K);A  = @(x) real(ifft2(K.*fft2(x)));      % Forward operatorAH = @(x) real(ifft2(KC.*fft2(x)));     % Inverse operatorgrad_f_x = @(x) AH(A(x) - Y);           % Gradient operatorxx = ones(size(Y));yy = 0*ones(size(Y));Lip = norm(grad_f_x(xx) - grad_f_x(yy), 2)/norm(xx - yy, 2); % A general %           calculation for Lipschitz constant via gradiend Lipschitz %           definition of data fidelity term.mu = 1.5/Lip; % error parameterPSNR_noisy = psnr(uint8(Y), uint8(X));RMSE_noisy = sqrt(mean((Y - X).^2));delta_x = inf;x_hat = zeros(size(Y));iter = 1;maxIter = 250;errorCriterion = 1e-3;old_X = x_hat;gamma = 17*sqrt(mu); % The lower bound for guaranteeing convergence sqrt(mu)/2tic;while (delta_x(iter) > errorCriterion) && (iter < maxIter)    iter = iter + 1;    Z = x_hat - (mu)*(grad_f_x(x_hat));    x_hat = CauchyProx(Z, gamma, mu);    delta_x(iter) = max(abs( x_hat(:) - old_X(:) )) / max(abs(old_X(:))); % Error calculation    old_X = x_hat;endtimeSim = toc;PSNR_regularized = psnr(uint8(x_hat), uint8(X));RMSE_regularized = sqrt(mean((x_hat(:) - X(:)).^2));figure;set(gcf, 'Position', [100 100 900 400])subplot('Position', [0.0101, 0.05001, 0.3, 0.95])imshow(uint8(X));       % Originaltitle('Original Image')subplot('Position', [0.3401, 0.05001, 0.3, 0.95])imshow(uint8(Y));       % Blurredtitle(['Blurred Image (PSNR = ' num2str(PSNR_noisy) ' dB)'])subplot('Position', [0.6701, 0.05001, 0.3, 0.95])imshow(uint8(x_hat));   % Cauchy Reconstructedtitle(['CPS Reconstructed (PSNR = ' num2str(PSNR_regularized) ' dB)'])fprintf('Cauchy proximal splitting (CPS) for image deblurring\nSolved after %d iterations in %.3f seconds\nNoisy PSNR = %.3f dB\nReconstructed PSNR = %.3f dB\n', iter, timeSim, PSNR_noisy, PSNR_regularized)

3 仿真结果

4 参考文献

[1] Karakus O ,  Mayo P ,  Achim A . Convergence Guarantees for Non-Convex Optimisation with Cauchy-Based Penalties[J].  2020.

博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。

​