Nonlinear diffusion by PDE in image processing

$$\frac{\partial u}{\partial t}=div(c(|\nabla u{|}^2)\nabla u)$$当$c=1$时是热方程。

Images are piecewise smooth functions. We want to smooth the image while preserving the sharp edges. 故取$$c(s)=\frac{1}{\sqrt{1+s}}$$

pic = imread('pic0.jpg');
u = imresize(pic,0.2);
u = rgb2gray(u);
imwrite(u,'pic.jpg');

pic = imread('pic.jpg');
u = double(pic);
s = size(u); M = s(1); N = s(2);
n_loop = 100;

% circular boundary value 
for x = 2:M-1
	u(x,1) = u(x,N-1);
	u(x,N) = u(x,1+1);
end
for y = 2:N-1
	u(1,y) = u(M-1,y);
	u(M,y) = u(1+1,y);
end

div_x = zeros(M,N); div_y = zeros(M,N);
for i = 1:n_loop
	for x = 2:M-1
		for y = 2:N-1
			div_x(x,y) = (u(x,y)-u(x-1,y))/((1+(u(x,y)-u(x-1,y))^2+(u(x,y)-u(x,y-1))^2)^0.5);
			div_y(x,y) = (u(x,y)-u(x,y-1))/((1+(u(x,y)-u(x-1,y))^2+(u(x,y)-u(x,y-1))^2)^0.5);
		end
	end
	for x = 3:M-1
		for y = 3:N-1
			u(x,y) = u(x,y)+div_x(x,y)-div_x(x-1,y)+div_y(x,y)-div_y(x,y-1);
		end
	end
	disp(i);
	if mod(i,10) == 0
		string = ['u_',num2str(i),'.jpg'];
		imwrite(uint8(u),string);
	end
end