基于偏微分方程去噪-热传导模型

来源：互联网发布：淘宝网卖假货会坐牢吗编辑：程序博客网时间：2024/05/09 15:50

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

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

0 0