1热传导方程
假设图像属于有界变差空间，那么，噪声图像应该满足两个条件：（1）噪声图像和原始图像相差不是特别大；（2）原始图像属于有界变差空间，那么，通过图像去噪可以建立为求解如下能量泛函的最优解的问题：
该能量泛函对应的Euler-Lagrange方程如下：
利用最速下降法，上述Euler-Lagrange方程可以转化为如下的PDEs的初边值问题：
该算法实现比较简单，matlab代码如下，但是，该模型实际上时对图像进行的模糊（注意到其中的laplace算子）
clear all;
close all;
clc;
Io=imread('picture.jpg'); %读入图像
if(ndims(Io)==3) %维数
Io=rgb2gray(Io); %转成灰色 rgb表示红绿蓝
end
Io=double(Io);
std_n=10; %高斯噪声标准差
var_n=std_n^2; %高斯噪声标准差
NI=randn(size(Io))*std_n; In=Io+NI; %在原图像上加噪声
dt=0.25; %网比（一般对于n维
%dt<= (1/2)^n这样子差分方程
%迭代才稳定）
N=100; %迭代次数
lambda=0; %lambda赋初值
[Max_J1 Max_J2 Min_J3 ALLPSNR ALLSNR ALLMAE J]=HeatEq(In,Io,dt,N,lambda); %方程迭代（热方程迭代）
[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]=HeatEq(In,Io,dt,N,lambda)
% 其中Io表示无噪声的原始图像
% In表示带高斯噪声的图像
% dt表示PED差分实现时的时间步长，一般较小(例如0.25，步长与差分方程的稳定性相关，这里不详细解释）
% N表示迭代次数
% lambda表示权重系数，越大表示去噪的图像和噪声图像相差应该越小
J = In;
Max_J1 = J;
Max_J2 = J;
Min_J3 = J;
for i=1:N
i
DxxI=J([2:end end],:)+J([1 1:end-1],:)-2*J; %函数关于X方向求二阶偏导
DyyI=J(:,[2:end end])+J(:,[1 1:end-1])-2*J; %函数关于Y方向求二阶偏导
J=J+dt*(DxxI+DyyI)-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
