如何计算被柯西噪声破坏的去噪图像的psnr值

Question

下面是我对一篇去除柯西噪声的论文的实现。噪声图像的psnr值为19，但当我计算恢复图像un的psnr值时，它返回17，这甚至小于噪声图像的psnr值。我想我的代码在psnr计算中有问题。

%%
clear memory;
clear all
close all
clc;
%% Initialization
refimg = im2double(imread('cameraman256.png')); % original image
img_height = size(refimg,1);
img_width = size(refimg,2);
refimg = refimg(1:img_height,1:img_width);
padNum = 5;
refimg = padarray(refimg,[padNum,padNum],'symmetric');
[mm,nn]=size(refimg);
img_height = size(refimg,1);
img_width = size(refimg,2);
%% Producing the degraded image
A = 1;  % A =1 for image denoising
sz = size(refimg);
rng(0);
r1 = randn(sz); % (using randn because I don't have the statistics toolbox)
r2 = randn(sz);
n = 0.02; % the noise level
u0 = refimg + n.*(r1./r2);
u0 = min(u0,1); % clamp large values to 1
u0 = max(u0,0); % clamp small values to 0
figure(1); imshow(u0(padNum+1:mm-padNum,padNum+1:nn-padNum),'border','tight');
%% Initial values for the primal-dual algorithm
tol = 1e-3;
nIter = 1e3;
options.order = 1; options.bound = 'sym';
un = u0;
wn = u0;
bun = un;
bwn = wn;
pxn = zeros(ny,nx);
pyn = zeros(ny,nx);
q = zeros(ny,nx);
[gxn,gyn] = grad(u0,options);
bgxn = gxn;
bgyn = gyn;
gamma = sqrt(2)/10;
lambda = 0.7;
mu = 6.25; 
tau = 0.3;
sigma = 0.3;
%% Primal-dual Algorithm
for j = 1:nIter 
%%%%%%%%%%solve the subproblem p
[ux,uy]=grad(bun,options);
pxn = pxn+sigma*(bgxn-ux);
pyn = pyn+sigma*(bgyn-uy);

%%%%%%%%%%%solve the subproblem q
AUk = Au(bun);
q = q+sigma*(bwn-AUk);
%%%%%%%%%%solve the subproblem g
goldxn = gxn;
goldyn = gyn;
txn = gxn-tau*pxn;
tyn = gyn-tau*pyn;
sn = max(1e-6,sqrt(txn.^2+tyn.^2));
gxn = txn./sn.*max(0,sn-tau);
gyn = tyn./sn.*max(0,sn-tau);

%%%%%%%%%%%solve the subproblem wn
u_medfilter = medfilt2(u0);
wold = wn;

a = mu.*lambda.*tau+1;
b = -(mu.*lambda.*tau.*(2.*u0+u_medfilter)-tau.*q+2.*u0+wold);
c = tau.*lambda+mu.*lambda.*tau.*(gamma.^2+u0.^2+2.*u_medfilter.*u0)-2.*tau.*q.*u0+...
gamma.^2+u0.^2+2.*wold.*u0;
d = -tau.*lambda.*u0-mu.*lambda.*tau.*u_medfilter.*(gamma.^2+u0.^2)+tau.*q.*(gamma.^2+u0.^2)...
-wold.*(gamma.^2+u0.^2);

qval = (3.*a.*c-(b.^2))./(9.*(a.^2));
rval = (9.*a.*b.*c-27.*(a.^2).*d-2.*(b.^3))./(54.*(a.^3));
deltaval = qval.^3+rval.^2;

wn = nthroot(rval+real(sqrt(deltaval)),3)+nthroot(rval-real(sqrt(deltaval)),3)-((b)./(3.*a));

%%%%%%%%%%%solve the subproblem un
uold = un;
Asqk = Atu(q);
un = un+tau*(Asqk-div(pxn,pyn,options)); % the restored image

%      un = min(1,max(0.01,un));

if (norm(un-uold, 'fro')/norm(uold,'fro')<tol)  
break;
end 

bun = 2*un-uold;
bwn = 2*wn-wold;
bgxn = 2*gxn-goldxn;
bgyn = 2*gyn-goldyn;

%      PSNR_restoredimage = psnr(refimg(padNum+1:mm-padNum,padNum+1:nn-b 
padNum),un(padNum+1:mm-padNum,padNum+1:nn-padNum))
end
%% Dispaly results
New_un = un;
refimg = refimg(padNum+1:mm-padNum,padNum+1:nn-padNum);
u0 = u0(padNum+1:mm-padNum,padNum+1:nn-padNum);
u_medfilter = u_medfilter(padNum+1:mm-padNum,padNum+1:nn-padNum); 
New_un = New_un(padNum+1:mm-padNum,padNum+1:nn-padNum);   
PSNR_noisy = psnr(refimg,u0)
PSNR_med = psnr(refimg,u_medfilter)
PSNR_restoredimage = psnr(New_un,refimg)
figure(2); imshow([refimg,u_medfilter,New_un],'border','tight'); 
%%

Answer 1

问题不在psnr计算中。在这个问题中，我的代码的问题是，我正在考虑去噪的情况，但我忘记了在原对偶算法中省略模糊核。此外，与去模糊情况相对应的初始参数必须为零。以下是已更正的代码。现在psnr值为28，与论文类似。

%%
clear;
%% Initialization
refimg = im2double(imread('cameraman256.png')); % original image
img_height = size(refimg,1);
img_width = size(refimg,2);
refimg = refimg(1:img_height,1:img_width);
%% Producing the noisy image
sz = size(refimg);
rng(0);
r1 = randn(sz); % (using randn because I don't have the statistics toolbox)
r2 = randn(sz);
n = 0.02; % the noise level
u0 = refimg + n.*(r1./r2);
u0 = min(u0,1); % clamp large values to 1
u0 = max(u0,0); % clamp small values to 0
figure(1); imshow(u0,'border','tight');
%% Initial values for the primal-dual algorithm
init = u0;
[ny,nx] = size(init);
tol = 1e-3;
nIter = 1e3;
options.order = 1; options.bound = 'sym';
un = u0;
wn = zeros(ny,nx);
bun = un;
bwn = zeros(ny,nx);
pxn = zeros(ny,nx);
pyn = zeros(ny,nx);
q = zeros(ny,nx);
[gxn,gyn] = grad(u0,options);
bgxn = gxn;
bgyn = gyn;
gamma = sqrt(2)/10;
lambda = 0.7;
mu = 6.25; 
tau = 0.3;
sigma = 0.3;
%% Primal-dual Algorithm
for j = 1:nIter 
%%%%%%%%%%solve the subproblem p
[ux,uy]=grad(bun,options);
pxn = pxn+sigma*(bgxn-ux);
pyn = pyn+sigma*(bgyn-uy);

%%%%%%%%%%%solve the subproblem q
AUk = bun;
q = q+sigma*(bwn-AUk);
%%%%%%%%%%solve the subproblem g
goldxn = gxn;
goldyn = gyn;
txn = gxn-tau*pxn;
tyn = gyn-tau*pyn;
sn = max(1e-6,sqrt(txn.^2+tyn.^2));
gxn = txn./sn.*max(0,sn-tau);
gyn = tyn./sn.*max(0,sn-tau);

%%%%%%%%%%%solve the subproblem wn
u_medfilter = medfilt2(u0);
wold = wn;

a = mu.*lambda.*tau+1;
b = -(mu.*lambda.*tau.*(2.*u0+u_medfilter)-tau.*q+2.*u0+wold);
c = tau.*lambda+mu.*lambda.*tau.*(gamma.^2+u0.^2+2.*u_medfilter.*u0)-2.*tau.*q.*u0+...
gamma.^2+u0.^2+2.*wold.*u0;
d = -tau.*lambda.*u0-mu.*lambda.*tau.*u_medfilter.*(gamma.^2+u0.^2)+tau.*q.*(gamma.^2+u0.^2)...
-wold.*(gamma.^2+u0.^2);

qval = (3.*a.*c-(b.^2))./(9.*(a.^2));
rval = (9.*a.*b.*c-27.*(a.^2).*d-2.*(b.^3))./(54.*(a.^3));
deltaval = qval.^3+rval.^2;

wn = nthroot(rval+real(sqrt(deltaval)),3)+nthroot(rval-real(sqrt(deltaval)),3)-((b)./(3.*a));

%%%%%%%%%%%solve the subproblem un
uold = un;
Asqk = q;
un = un+tau*(Asqk-div(pxn,pyn,options)); % the restored image


if (norm(un-uold, 'fro')/norm(uold,'fro')<tol)  
break;
end 

bun = 2*un-uold;
bwn = 2*wn-wold;
bgxn = 2*gxn-goldxn;
bgyn = 2*gyn-goldyn;

end
%% Dispaly results
PSNR_noisy = psnr(u0,refimg)
PSNR_med = psnr(u_medfilter,refimg)
PSNR_restoredimage = psnr(un,refimg)
figure(2); imshow([refimg,u_medfilter,un],'border','tight'); 
%%