%求解四階偏微分方程
%%%%%%%%%%%%%%%%%%
%準備工作：
%1.時間tao和空間h的劃分
clear all
h=1.2;
tao=0.2;
h2=1.0/(h*h);
n=0;%記錄循環次數
p=3;q=3;%矩陣的大小
%2.差分格式的乘積矩陣
A=[0 1 0;1 -4 1;0 1 0];
B=zeros(3);
C=[1 1 1;1 1 1;1 1 1];
%%%%%%%%%%%%%%%%%%%%%
% 讀入圖像
[filename,filepath]=uigetfile('*.jpg','請輸入一個圖像');
img=strcat(filepath,filename);

%圖像灰度處理
color=imfinfo(img);
colortype=color.ColorType;   %獲取圖象顏色類型
%類型判斷和相應的讀取圖像
switch(colortype)
    case 'truecolor'
        im0=rgb2gray(imread(img));%轉化為灰度圖像
    case  'indexed'
        [im0,map]=imread(img);
        im0=ind2gray(im,map);     %轉化為灰度圖像
    otherwise
        im0=imread(img);
end
%%%%%%%%%%%%%%%%%%%%%%
% gaussian
% im0 = imnoise(im0,'gaussian');
figure(1),imshow(im0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%開辟所需空間,把原圖像數字矩陣擴充兩行兩列
%1.準備工作
% im0=rgb2gray(im0);
[m,n]=size(im0);%求取矩陣大小
im1=zeros(m+2,n+2);%開辟一個新的矩陣空間
im2=zeros(m+2,n+2);%開辟一個新的矩陣空間
im3=zeros(m+2,n+2);%開辟一個新的矩陣空間
mean0=zeros(m+2,n+2);%開辟一個新的矩陣空間
var0=zeros(m+2,n+2);%開辟一個新的矩陣空間
max0=zeros(m+2,n+2);%開辟一個新的矩陣空間
g0=zeros(m+2,n+2);%開辟一個新的矩陣空間
g1=zeros(m+2,n+2);%開辟一個新的矩陣空間
g2=zeros(m+2,n+2);%開辟一個新的矩陣空間
%2.實現
im1(2:m+1,2:n+1)=im0(:,:);%把原來的矩陣值賦給它
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while 1
    %求矩陣g1
    %1.求im0的5五點差分格式，得到矩陣im1
    for i=2:m+1
        for j=2:n+1
            %
            B=im1(i-1:i+1,j-1:j+1);
            im2(i,j)=h2*sum(sum(A.*B));
            %計算均值
            mean0(i,j)=sum(sum(B))/(p*q);%3*3
            %計算方差
            B1=B-mean0(i,j);
            var0(i,j)=sum(sum(B1.^2))/(p*q);
        end
    end
    %2.求g0=cu
        max0=max(max(var0));
        g0=1./(var0./max0.*abs(im2)/255+1);
    %3.求g0的5五點差分格式，得到矩陣g1
    for i=2:m+1
        for j=2:n+1
            B=g0(i-1:i+1,j-1:j+1);
            g1(i,j)=h2*sum(sum(A.*B));
        end
    end
    %4.求最后的差分格式
    g2=tao*g1;
    im3=im1-g2;%important
    %滿足條件就退出循環
    n=n+1;
    if n>10
        break;
    end
    %更新原始矩陣
    im1=im3;
end %while
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%輸出最后的矩陣圖像
figure(2),imshow(mat2gray(im3(2:m,2:n)));