tic;
f=imread('rice.tif');  %灰度圖
% f=double(a);           %卷積運算不支持uint8類型？？
TPh=zeros(1,256); %非零梯度象素的總數
TPl=zeros(1,256);
EAGhigh=zeros(1,256);%圖像中非零梯度象素的平均梯度(有效平均梯度)
EAGlow=zeros(1,256);
h1=[-1 0 1;-2 0 2;-1 0 1];
h2=[1 2 1;0 0 0;-1 -2 -1];
fhigh=f;
flow=f;
for x=0:255;
 fhigh(f>=x)=x;        %高端剪切的梯度圖
 flow(f<=x)=x;         %低端剪切的梯度圖
 gh1=conv2(fhigh,h1,'valid');
 gh2=conv2(fhigh,h2,'valid');%不考慮邊界補零，返回矩陣小于原矩陣
 gl1=conv2(flow,h1,'valid');
 gl2=conv2(flow,h2,'valid');
%  gh=sqrt((gh1.*gh1)+(gh2.*gh2));
%  gl=sqrt((gl1.*gl1)+(gl2.*gl2));
 gh=abs(gh1)+abs(gh2);   
 gl=abs(gl1)+abs(gl2);
 TGh(x+1)=sum(sum(gh));  %梯度圖的總梯度值
 TGl(x+1)=sum(sum(gl));
 ph=(gh>0);
 pl=(gl>0);
 TPh(x+1)=sum(sum(ph));   %非零梯度象素的總數
 TPl(x+1)=sum(sum(pl));
 if  TPh(x+1)==0
        EAGhigh(x+1)=0;
    else
        EAGhigh(x+1)=TGh(x+1)/TPh(x+1);
    end
    if  TPl(x+1)==0;
        EAGlow(x+1)=0;
    else
        EAGlow(x+1)=TGl(x+1)/TPl(x+1);    %有效平均梯度
    end       
end
[EAGh,Lhigh]=max(EAGhigh);
[EAGl,Llow]=max(EAGlow);
Lhigh,Llow,
 toc  %19.478s