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
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的的方差
if var_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));

if nostotal == 0
    s = 999.99; %% INF. clean image
else
    s = 10 * log10(peak./nostotal);%计算峰值性噪比
end
return

转载地址：https://blog.csdn.net/hit1524468/article/details/53444812?locationNum=12&fps=1