对应的Euler-Lagrange方程同前所述一样，如下：
由最速下降法获得上述能量泛函对应的PDEs初边值问题如下：
下面为该偏微分方程的matlab实现代码：
clear all
close all
clc
Io=imread('picture.jpg');% 读入一幅图像
% 选择一个颜色矩阵,并且变成双浮点型的,这一步没有可能报错.
if(ndims(Io) == 3)
Io = rgb2gray(Io);
end;
Io=double(Io);
% % % %%% Add noise %%% % % %
std_n=20; % 高斯噪声比准差
var_n=std_n^2; % 高斯噪声比准差
NI = randn(size(Io))*std_n; % 白色高斯噪声
In = Io + NI; % 把噪声加到原图上面
N =100; %迭代图像
dt=0.2; %网比（一般对于n维时，dt<= (1/2)^n这样子差分方程迭代才稳定）
lambda=0.01; %给lambda赋初值
tic
[Max_J1 Max_J2 Min_J3 ALLPSNR ALLSNR ALLMAE J] =TV(In,Io,dt,N,lambda,var_n);%调用函数
toc
[MaxPSNR, Index1]=max(ALLPSNR)
[MaxSNR, Index2]=max(ALLSNR)
[MinMAE, Index3]=min(ALLMAE)
Print(Io,In,J,Max_J1,Max_J2,Min_J3,ALLPSNR,ALLSNR,ALLMAE,N);
function [Max_J1 Max_J2 Min_J3 ALLPSNR ALLSNR ALLMAE J]=TV(In,Io,dt,N,lambda,var_n)
J=In;
Max_J1 = J;
Max_J2 = J;
Min_J3 = J;
ep=0.0001;
for i=1:N
i
DfJx=J([2:end end],:)-J; %函数关于X的一阶偏导（向后差分）
DbJx=J-J([1 1:end-1],:); %函数关于X的一阶偏导（向前差分）
DfJy=J(:,[2:end end])-J; %函数关于Y的一阶偏导（向后差分）
DbJy=J-J(:,[1 1:end-1]); %函数关于Y的一阶偏导（向前差分）
TempDJx=(ep+DfJx.*DfJx+((sign(DfJy)+sign(DbJy)).*min(abs(DfJy),abs(DbJy))./2).^2).^(1/2);%求梯度的模
TempDJy=(ep+DfJy.*DfJy+((sign(DfJx)+sign(DbJx)).*min(abs(DfJx),abs(DbJx))./2).^2).^(1/2);
DivJx=DfJx./TempDJx;
DivJy=DfJy./TempDJy;
% DivJx=DfJx./TempDJx^p;
% DivJy=DfJy./TempDJy^p;
%求散度
Div=DivJx-DivJx([1 1:end-1],:)+DivJy-DivJy(:,[1 1:end-1]);
% update lambda (fidelity term)
lambda = max(mean(mean(Div.*(J-In)))./var_n,0)
J= J+ dt * Div -dt*lambda*(J-In); %产生迭代
NowPSNR = psnr(uint8(J),Io) %调用psnr函数
ALLPSNR(i)=NowPSNR;
if i>1 && ALLPSNR(i-1) < ALLPSNR(i)
Max_J1 = J;
end
NowSNR = snr(uint8(J),Io)
ALLSNR(i) = NowSNR;
if i>1 && ALLSNR(i-1) < ALLSNR(i)
Max_J2 = J;
end
NowMAE = mae(uint8(J),Io)
ALLMAE(i) = NowMAE;
if i>1 && ALLMAE(i-1) > ALLMAE(i)
Min_J3 = J;
end
end
function Print(Io,In,J,Max_J1,Max_J2,Min_J3,ALLPSNR,ALLSNR,ALLMAE,N)
figure(1)
subplot(2,2,1)
imshow(Io,[]);
title('原图像')
subplot(2,2,2)
imshow(In,[]);
title('加噪声之后的图像')
subplot(2,2,3)
imshow(Io,[]);
title('原图像')
subplot(2,2,4)
imshow(J,[]);
title('处理结果')
figure(2)
subplot(2,2,1)
imshow(Max_J1,[]);
title('ALLPSNR值最大时图像')
subplot(2,2,2)
imshow(Max_J2,[]);
title('ALLSNR值最大时图像')
subplot(2,2,3)
imshow(Min_J3,[]);
title('ALLPMAE值最小时图像')
subplot(2,2,4)
imshow(J,[]);
title('处理结果')
x=1:N;
figure(3)
subplot(2,2,1)
plot(x,ALLPSNR)
title('ALLPSNR图像')
subplot(2,2,2)
plot(x,ALLSNR)
title('ALLSNR图像')
subplot(2,2,3)
plot(x,ALLMAE)
title('ALLPMAE图像')
[Ny,Nx]=size(J);
x=1:Nx;
level=fix(Ny/2);
y=J(level,:);
y1=Io(level,:);
y2=In(level,:);
figure(4)
subplot(2,1,1); plot(x,y,x,y1);
title('SmoothImage And OriginalImage')
subplot(2,1,2); plot(x,y,x,y1,x,y2);
title('NoiseImage And OriginalImage')
function s = snr(noisydata, original)
%将noisydata,original转化为double型
noisydata = double(noisydata);
original = double(original);
mean_original = mean(original(:));%求original的平均值
tmp = original - mean_original;
var_original = sum(sum(tmp.*tmp));%求original的方差
noise = noisydata - original;%求noise的平均值
mean_noise = mean(noise(:));
tmp = noise - mean_noise;
var_noise = sum(sum(tmp.*tmp));%求noise的的方差
ifvar_noise == 0
s = 999.99; %% INF. clean image
else
s = 10 * log10(var_original / var_noise);%compute signal-to-noise-ratio (SNR) of a noisy signal/image
end
return
function E = mae(noisydata, original)
%将noisydata,original转化为double型
noisydata=double(noisydata);
original=double(original);
[m,n] = size(noisydata);
noise = abs(noisydata - original);
nostotal = sum(sum(noise));
E=nostotal/(m*n);%compute root-mean-square-error (RMSE) of a noisy signal/image
Return
function s = psnr(noisydata, original)
%将noisydata,original转化为double型
noisydata=double(noisydata);
original=double(original);
[m,n] = size(noisydata);%获得noisydata矩阵的行数与列数
peak=255*255*m*n;%计算峰值
noise =noisydata - original;
nostotal = sum(sum(noise.*noise));
ifnostotal == 0
s = 999.99; %% INF. clean image
else
s = 10 * log10(peak./nostotal);%计算峰值性噪比
end
return