fangli.m

GM(1,1)模型 MATLAB

%修正p值gm(1,1) 检测
clc;clear;
m2=imread('lena_gray.jpg');
m=double(m2)./255;
mm2=m;
p=0.5;
[n1,n2]=size(m);
for i=1:n1-1                            %水平方向
    for j=1:n2-2
        AA2=[m(i,j);m(i,j)+m(i,j+1);m(i,j)+m(i,j+1)+m(i,j+2);m(i,j)+m(i,j+1)+m(i,j+2)+m(i+1,j)];
        %A2=[AA2(1),AA2(1)+AA2(2),AA2(1)+AA2(2)+AA2(3),AA2(1)+AA2(2)+AA2(3)+AA2(4)];
        B2=[-0.5*(AA2(1)+AA2(2)),1;-0.5*(AA2(2)+AA2(3)),1;-0.5*(AA2(3)+AA2(4)),1];
        F2=[m(i,j+1);m(i,j+2);m(i+1,j)];	%AA2为累加过的，这里要不累加的

        E2=inv(B2'*B2)*B2'*F2;      %E2(1)会变得非常小

       % G2=(m(i,j)-E2(2)/E2(1))*exp(-4*E2(1)/p)+E2(2)/E2(1);
       % T2=(m(i,j)-E2(2)/E2(1))*exp(-3*E2(1)/p)+E2(2)/E2(1);
       % S2=(m(i,j)-E2(2)/E2(1))*exp(-2*E2(1)/p)+E2(2)/E2(1);
       % mm2(i+1,j+1)=G2-2*T2+S2;
		mm2(i+1,j+1) = ( (1-exp(E2(1)*p)) * (m(i,j) - E2(2)/E2(1)) * exp(-E2(1)*4/p) ); %此为化简后的式子，减少不必要的计算
	   %		R2=uint8(round(mm2*255));       %可以提取到循环外部
    end
end
R2=uint8(round(mm2*255)); 
imwrite(R2,'xiesong2.jpg');figure;imshow(R2);