% AMSS nonlinear diffusion

function y=AMSS(x,n)

x=padarray_new(x,1);
[M,N]=size(x);
IRAC8 = 0.35355339;
IRAC2P2 = 0.29289322 ;
IRAC2 =   0.70710678 ;
RAC8P4 = 6.8284271;
step=0.125;
y=x;

for iter=1:n
    for i=2:M-1
        for j=2:N-1
            c1=x(i+1,j+1)-x(i-1,j-1);
            d1=x(i-1,j+1)-x(i+1,j-1);
            ax = x(i+1,j)-x(i-1,j)+IRAC8*(c1-d1); 
            ay = x(i,j+1)-x(i,j-1)+IRAC8*(c1+d1);
            axy = ax*ay;
      ax = ax*ax;
      ay = ay *ay;
      an2 = ax+ay;
      if an2>0
      ax =ax / an2; ay=ay / an2; axy=axy / an2;
	l0 = 0.5 - axy*axy;	
	GK =     x(i,j)    * ( -4.0*l0 ) ...
	    + (x(i-1,j)+x(i+1,j)) * ( 2.0*l0 -      ax              ) ...
	    + (x(i,j-1)+x(i,j+1)) * ( 2.0*l0 -           ay         ) ...
	    + (x(i+1,j+1)+x(i-1,j-1)) * ( -l0 + 0.5 * ( ax + ay - axy ) ) ...
	    + (x(i+1,j-1)+x(i-1,j+1)) * ( -l0 + 0.5 * ( ax + ay + axy ) );
    GK = GK* an2*IRAC2P2*IRAC2P2;
    % if power==1
	  if (GK>=0.0) 
	    % here we use the fact that |Du|K^(1/3) = (|Du|^2*GK)^(1/3) 
	    y(i,j) = x(i,j) + step * (GK^(1/3));
	  else 
	    y(i,j) = x(i,j) - step * ((-GK)^(1/3));
        end
    else
    y(i,j)=x(i,j);
    end
end
    
end
x=y;
end
y=padarray_new(y,-1);
        % else 
	  % MCM diffusion 
	  % y(i,j) = x(i,j) + step * GK;