function MCM(image_name,t,out_filname)
% AMSS
% image_name  原圖像文件名
% t           時間
% out_filname 輸出文件名

if nargin<2
    disp('輸入參數錯誤!')
    return;
end
if t<0.1
    t=0.1;
end
if isempty(out_filname)
    out_filname=strcat('MCM_',image_name);
end
if nargin<2
    disp('輸入參數錯誤!')
    return;
end
if t<0.1
    t=0.1;
end

constant=0.35355339;

image_I=imread(image_name);
subplot(1,2,1);
imshow(image_I,[]);
[size_r,size_c]=size(image_I);

   C1=zeros(size_r,size_c);
   D1=zeros(size_r,size_c);

   U_dx=zeros(size_r,size_c);
   U_dy=zeros(size_r,size_c);
   Gk=zeros(size_r,size_c);
 
 for inde=1:t/0.1
    image_ex=matrix_expand(image_I);
    
    for row=1:size_r
        for col=1:size_c 
            C1(row,col)=image_ex(row+2,col+2)-image_ex(row,col);
            D1(row,col)=image_ex(row,col+2)-image_ex(row+2,col);
        end
    end
            
    for row=1:size_r
        for col=1:size_c  
            U_dx(row,col)=image_ex(row+2,col+1)-image_ex(row,col+1);
            U_dy(row,col)=image_ex(row+1,col+2)-image_ex(row+1,col);
        end
    end
    
    U_dx=U_dx+constant*(C1-D1);
    U_dy=U_dy+constant*(C1+D1);
    
    a_x=U_dx.*U_dx;
    a_y=U_dy.*U_dy;
    a_xy=U_dx.*U_dy;
    a_s2=a_x+a_y;
    
    for row=1:size_r
        for col=1:size_c 
            if a_s2(row,col)~=0
                a_x(row,col)=a_x(row,col)/a_s2(row,col);
                a_y(row,col)=a_y(row,col)/a_s2(row,col);
                a_xy(row,col)=a_xy(row,col)/a_s2(row,col);
            end
        end
    end
    
    L_0=0.5-a_xy.*a_xy;
    
    beta_0=-4*L_0;
    beta_1=2*L_0-a_x;
    beta_2=2*L_0-a_y;
    beta_3=-L_0+0.5*(a_x+a_y-a_xy);
    beta_4=-L_0+0.5*(a_x+a_y+a_xy);
    
    for row=1:size_r
        for col=1:size_c 
            Gk(row,col)=image_ex(row+1,col+1)*beta_0(row,col)+( image_ex(row,col+1)+image_ex(row+2,col+1) )*beta_1(row,col)+( image_ex(row+1,col)+image_ex(row+1,col+2) )*beta_2(row,col)+( image_ex(row,col)+image_ex(row+2,col+2) )*beta_3(row,col)+( image_ex(row,col+2)+image_ex(row+2,col))*beta_4(row,col);
        end
    end
    
    v=double(image_I)+0.1*Gk;
    image_I=uint8(v);
end

  subplot(1,2,2);
  imshow(uint8(v));  
  imwrite(uint8(v),out_filname);    
% －－－－－－－－－－－－－－－－－－－－－－end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function matrix_ex=matrix_expand(matrix)
% 對矩陣進行延拓   
   [size_r,size_c]=size(matrix);
   matrix_ex=zeros(size_r+2,size_c+2);
   % 中
   for row=1:size_r
       for col=1:size_c
           matrix_ex(row+1,col+1)=matrix(row,col);
       end
  end
  % 左，右
  for row=1:size_r
      matrix_ex(row+1,1)=matrix(row,1);
      matrix_ex(row+1,size_c+2)=matrix(row,size_c);
  end
  % 上，下
  for col=1:size_c
      matrix_ex(1,col+1)=matrix(1,col);
      matrix_ex(size_r+2,col+1)=matrix(size_r,col);
  end
  matrix_ex(1,1)=matrix(1,1);
  matrix_ex(1,size_c+2)=matrix(1,size_c);

  matrix_ex(size_r+2,1)=matrix(size_r,1);
  matrix_ex(size_r+2,size_c+2)=matrix(size_r,size_c);
% -------------------------------------------------