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