function?Jd=diffusion（JmethodNKdtsigma2）
%?private?function:?diffusion?（by?Guy?Gilboa）:
%?Jd=diffusion（JmethodNK）
%?Simulates?N?iterations?of?diffusion?parameters:
%?J?=??source?image?（2D?gray-level?matrix）?for?diffusio
%?method?=??‘lin‘:??Linear?diffusion?（constant?c=1）.
%???????????‘pm1‘:?perona-malik?c=exp{-（|grad（J）|/K）^2}?[PM90]
%???????????‘pm2‘:?perona-malik?c=1/{1+（|grad（J）|/K）^2}?[PM90]
%???????????‘rmp‘:?complex?valued?-?ramp?preserving?[GSZ01]
%?K????edge?threshold?parameter
%?N????number?of?iterations
%?dt???time?increment?（0?%?sigma2?-?if?present?calculates?gradients?of?diffusion?coefficient
%??????????convolved?with?gaussian?of?var?sigma2?（Catte?et?al?[CLMC92]）

if?~exist（‘N‘）
???N=1;
end
if?~exist（‘K‘）
???K=1;
end
if?~exist（‘dt‘）
???dt=0.2;
end
if?~exist（‘sigma2‘）
???sigma2=0;
end

[NyNx]=size（J）;?

if?（nargin<3）?
???error（‘not?enough?arguments?（at?least?3?should?be?given）‘）;
end

for?i=1:N;???
???%?gaussian?filter?with?kernel?5x5?（Catte?et?al）
???if?（sigma2>0）?
??????Jo?=?J;???%?save?J?original
??????J=gauss（J5sigma2）;??
???end

	%?calculate?gradient?in?all?directions?（NSEW）
	In=[J（1:）;?J（1:Ny-1:）]-J;
	Is=[J（2:Ny:）;?J（Ny:）]-J;
	Ie=[J（:2:Nx）?J（:Nx）]-J;
???Iw=[J（:1）?J（:1:Nx-1）]-J;

	%?calculate?diffusion?coefficients?in?all?directions?according?to?method
???if?（method?==?‘lin‘）
	???Cn=K;?Cs=K;?Ce=K;?Cw=K;
???elseif?（met