Effect of Gaussian Smoothing to mean image intensity

Question

I am implementing Perona-Malik Diffusion in using the formulas:

Creating image at time k+1 using equation

This is my implementation:

function gij=g(duij,lambda)
duij=-(duij.^2/lambda^2);
gij=exp(duij);
end

function deltauij=gradMag(u)
xx=[u(:,2:end),u(:,end)]-[u(:,1),u(:,1:end-1)];      
yy=[u(2:end,:);u(end,:)]-[u(1,:);u(1:end-1,:)];
xx=(xx.^2)/4;
yy=(yy.^2)/4;
deltauij = xx+yy;
deltauij = sqrt(deltauij);

function u=nld(u,delt,lambda,sigma)

filt= fspecial('gaussian',[3*ceil(sigma),3*ceil(sigma)],sigma);
utemp=imfilter(u,filt);
deltauij=gradMag(utemp);
gm=g(deltauij,lambda);
uxyplus=[u(:,2:end),u(:,end)];
uxyminus=[u(:,1),u(:,1:end-1)];
uxplusy=[u(2:end,:);u(end,:)];
uxminusy=[u(1,:);u(1:end-1,:)];

gxyplus=[gm(:,2:end),gm(:,end)];
gxyminus=[gm(:,1),gm(:,1:end-1)];
gxplusy=[gm(2:end,:);gm(end,:)];
gxminusy=[gm(1,:);gm(1:end-1,:)];

xx=0.5*(gxplusy+gm).*(uxplusy-u)-0.5*(gxminusy+gm).*(u-uxminusy);
yy=0.5*(gxyplus+gm).*(uxyplus-u)-0.5*(gxyminus+gm).*(u-uxyminus);
ddu=xx+yy;
del = delt*ddu;
u=min(255,max(0,del+u));

And I am calling iteratively I=nld_2(I,delta_t,lambda,sigma)
I am calculating mean intensity value of image in every iteration.The result is quite weird .Mean vs iteration plot is lik this

Since this is a smoothing operation I am expecting it to decrease mean intensity somehow but what I get is an increasing mean intensity.

Is this normal or I am missing anything?