用于图像恢复的即插即用 ADMM：定点收敛和应用附Matlab代码

原创于 2023-09-23 23:59:39 发布 · 363 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#matlab #开发语言

图像处理专栏收录该内容

871 篇文章

订阅专栏

本文介绍了一种具有定点收敛性的插件ADmm算法，用于解决图像恢复中的约束优化问题。研究了如何将不同去噪算法融入插件与播放Admm，探讨了保证收敛的条件，并展示了在超分辨率和单光子成像中的应用，对比了其与现有算法的性能。

✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知。

更多Matlab完整代码及仿真定制内容点击👇

智能优化算法神经网络预测雷达通信无线传感器电力系统

信号处理图像处理路径规划元胞自动机无人机

🔥 内容介绍

乘数交替法(ADmm)是解决图像恢复中约束优化问题的一种常用算法。在许多有用的特征中,Admm算法的一个关键特征是它的模块化结构,该结构允许我们插入任何现成的ADmm算法中子问题的图像去噪算法。由于插件的性质,这种类型的Admm算法被命名为"插件和播放Admm"。在最近的几篇论文中,"插件与播放"广告已经展示了有希望的实证结果。然而,目前还不清楚在什么条件下使用什么去噪算法来保证收敛性。另外,由于插件和播放Admm使用一种特定的方法来分解变量,所以对于常见的高斯和波斯图像恢复问题,是否可以快速实现还不清楚。
本文提出了一种具有可证明的定点收敛性的插用ADmm算法。我们证明,对于任何满足渐近准则的去噪算法,称为有界除数器,插入和播放的Admm收敛到一个固定点下的延续方案。我们还介绍了超分辨率和单光子成像的两个图像恢复问题的快速实现。我们比较了每一个问题类型中最先进的算法和插-玩算法,并给出了该算法的实验结果。

📣 部分代码

function out = PlugPlayADMM_general(y,A,lambda,method,opts)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%out = PlugPlayADMM_general(y,A,lambda,method,opts)%solves the general inverse problem%%inversion step: x=argmin_x(||Ax-y||^2+rho/2||x-(v-u)||^2)%denoising step: v=Denoise(x+u)%      update u: u=u+(x-v)%%Input:           y    -  the observed gray scale image%                 A    -  forward operator%              lambda  -  regularization parameter%              method  -  denoiser, e.g., 'BM3D'%       opts.rho       -  internal parameter of ADMM {1}%       opts.gamma     -  parameter for updating rho {1}%       opts.maxitr    -  maximum number of iterations for ADMM {20}%       opts.tol       -  tolerance level for residual {1e-4}   %       ** default values of opts are given in {}. %%Output:          out  -  recovered gray scale image %%Xiran Wang and Stanley Chan%Copyright 2016%Purdue University, West Lafayette, In, USA.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Check inputsif nargin<4    error('not enough input, try again \n');elseif nargin==4    opts = [];end% Check defaultsif ~isfield(opts,'rho')    opts.rho = 1;endif ~isfield(opts,'max_itr')    opts.max_itr = 20;endif ~isfield(opts,'tol')    opts.tol = 1e-4;endif ~isfield(opts,'gamma')    opts.gamma=1;endif ~isfield(opts,'print')    opts.print = false;end% set parametersmax_itr   = opts.max_itr;tol       = opts.tol;gamma     = opts.gamma;rho       = opts.rho;%initialize variablesdim         = size(y);          %dimension of the imageN           = dim(1)*dim(2);    %number of pixels in the imagev           = 0.5*ones(dim);x           = v;u           = zeros(dim);residual    = inf;%set function handle for denoiserswitch method    case 'BM3D'        denoise=@wrapper_BM3D;    case 'TV'        denoise=@wrapper_TV;    case 'NLM'        denoise=@wrapper_NLM;    case 'RF'        denoise=@wrapper_RF;    otherwise        error('unknown denoiser \n');end% main loopif opts.print==true    fprintf('Plug-and-Play ADMM --- General \n');    fprintf('Denoiser = %s \n\n', method);    fprintf('itr \t ||x-xold|| \t ||v-vold|| \t ||u-uold|| \n');enditr = 1;while(residual>tol&&itr<=max_itr)    %store x, v, u from previous iteration for psnr residual calculation    x_old = x;    v_old = v;    u_old = u;        G = @(z,trans_flag) gfun(z,trans_flag,A,rho,dim);    %inversion step    xtilde = v-u;    rhs    = [y(:); sqrt(rho)*xtilde(:)];    [x,~]  = lsqr(G,rhs,1e-3);    x      = reshape(x,dim);        %denoising step    vtilde = x+u;    vtilde = proj(vtilde);    sigma  = sqrt(lambda/rho);    v      = denoise(vtilde,sigma);        %update langrangian multiplier    u      = u + (x-v);        %update rho    rho = rho*gamma;        %calculate residual    residualx = (1/sqrt(N))*(sqrt(sum(sum((x-x_old).^2))));    residualv = (1/sqrt(N))*(sqrt(sum(sum((v-v_old).^2))));    residualu = (1/sqrt(N))*(sqrt(sum(sum((u-u_old).^2))));        residual = residualx + residualv + residualu;    if opts.print==true        fprintf('%3g \t %3.5e \t %3.5e \t %3.5e \n', itr, residualx, residualv, residualu);    end        itr = itr+1;endout = v;endfunction y = gfun(x,transp_flag,A,rho,dim)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% This gfun assumes that A is an N-by-N matrix% If A is M-by-N, then the following program can be changed% to x1 = x(1:M); x2 = (M+1:M+N)%% Stanley Chan% Nov 26, 2016%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%rows = dim(1);cols = dim(2);N    = rows*cols;if strcmp(transp_flag,'transp')         % y = A'*x    x1   = x(1:N);    x2   = x(N+1:2*N);    Atx  = A(x1,'transp');    y    = Atx + sqrt(rho)*x2;elseif strcmp(transp_flag,'notransp')   % y = A*x    Ax   = A(x,'notransp');    y    = [Ax; sqrt(rho)*x];endend