一、简介
1 用途
图像降噪(图像修复，复原中)，TVloss是一种较为有效的正则项，来保持图像的光滑性。
2 效果
二、源代码
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I=imread('toys.bmp'); % load image
I=double(I(20:120,10:105)); % cut a piece, convert to double
%%% Parameters
std_n=10; var_n=std_n^2; % Gaussian noise standard deviation
reduced_pw = 1.5*var_n; % power to reduce in first phase
sig_w = 5; ws=4*sig_w+1; % window size
%%%%%%%%%%%%%%
%%% Add noise
In = randn(size(I))*std_n;
I0 = I + In; % noisy input image
% show original and noisy images
figure(1); imshow(uint8(I)); title('Original')
figure(2); imshow(uint8(I0)); title('Noisy image')
snr_noisy=db(I,I0)
% run normal tv - strong denoising
tic;
J=I0;
ep_J=0.1; % minimum mean change in image J
J_old=0;
lam=0; iter=10; dt=0.2; ep=1;
while (mean(mean(abs(J - J_old))) > ep_J), % iterate until convergence
J_old=J;
J=tv(J,iter,dt,ep,lam,I0); % scalar lam
lam = calc_lam(J,I0,reduced_pw); % update lambda (fidelity term)
end % while
% figure(3); imshow(uint8(J)); title('residue TV')
% snr_residue= db(I,J)
Ir=I0-J; % Ir scalar
Pr = mean(mean(Ir.^2)); % power of residue
LV = loc_var(Ir,ws,sig_w^2); % local variance (local power of the residue )
% P=mean(mean(LV));
Pxy=1*(var_n^2)./LV; %% Sxy inverse proportional to the LV
%%% Varying Lambda
lamxy=zeros(size(I0));
J=I0; J_old=0;
ep_J=0.001;
%eps=0.01;
while (mean(mean(abs(J - J_old))) > ep_J), % iterate until convergence
J_old=J;
J=tv(J,iter,dt,ep,lamxy,I0); % adaptive lam
%J=tv(J,iter,dt,ep_J,lamxy,I0);
lamxy = calc_lamxy(J,I0,Pxy,sig_w); % update lambda (fidelity term)
end % while
figure(3); imshow(uint8(J)); title('Adaptive TV')
snr_adap= db(I,J)
toc
t1=toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Run scalar TV denoising for comparision
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
J=I0;
% params
ep_J = 0.01; % minimum mean change in image J
lam=0; J_old=0;
i=0;
while (mean(mean(abs(J - J_old))) > ep_J), % iterate until convergence
J_old = J;
J=tv(J,iter,dt,ep,lam,I0); % scalar lam
lam = calc_lam(J,I0,var_n,ep); % update lambda (fidelity term)
end % for i
% Ir=I0-J; % Ir scalar
% Pr = mean(mean(Ir.^2)); % power of residue
function Ig=gauss(I,ks,sigma2)
%private function: gauss (by Guy Gilboa):
% Ig=gauss(I,ks,sigma2)
% ks - kernel size (odd number)
% sigma2 - variance of Gaussian
[Ny,Nx]=size(I);
hks=(ks-1)/2; % half kernel size
if (Ny<ks) % 1d convolutin
x=(-hks:hks); %x:1x ks
flt=exp(-(x.^2)/(2*sigma2)); % 1D gaussian
flt=flt/sum(sum(flt)); % normalize
% expand
x0=mean(I(:,1:hks)); xn=mean(I(:,Nx-hks+1:Nx));%x0,xn :1 x hks;
eI=[x0*ones(Ny,ks) I xn*ones(Ny,ks)]; % ???
Ig=conv(eI,flt);
Ig=Ig(:,ks+hks+1:Nx+ks+hks); % truncate tails of convolution
else
%% 2-d convolution
x=ones(ks,1)*(-hks:hks); y=x';
flt=exp(-(x.^2+y.^2)/(2*sigma2)); % 2D gaussian
flt=flt/sum(sum(flt)); % normalize
% expand
if (hks>1)
xL=mean(I(:,1:hks)')'; xR=mean(I(:,Nx-hks+1:Nx)')';% xL,xR :Ny x 1
else
xL=I(:,1); xR=I(:,Nx);
end
eI=[xL*ones(1,hks) I xR*ones(1,hks)]; % Ny x Nx+2hks
if (hks>1)
xU=mean(eI(1:hks,:)); xD=mean(eI(Ny-hks+1:Ny,:)); % xU,xD: 1 x Nx+2hks
else
xU=eI(1,:); xD=eI(Ny,:);
end
三、运行结果
四、备注
版本：2014a
