pic2=imread（‘fig3.tif‘）;
subplot（231）;
imshow（pic2）;
title（‘原始圖像‘）;
[ROWCOL]=size（pic2）;???????????????????%取表示圖像的矩陣行列值數(shù)
count=imhist（pic2）;?????????????????????%直方圖統(tǒng)計????????????????????%圖像矩陣大小
N=ROW*COL;?????????????????????????????%圖像像素個數(shù)
L=256;?????????????????????????????????%指定圖像灰度級為256級
count=count/N;?????????????????????????%各灰度級出現(xiàn)的概率
%================尋找出現(xiàn)概率不為0的最小灰度值=============
for?i=2:L??????????????????????????
????if?count（i）~=0
????????st=i-1;
????????break;
????end
end
%================尋找出現(xiàn)概率不為0的最大灰度值=============
for?i=L:-1:1
????if?count（i）~=0
????????nd=i-1;
????????break;
????end
end
%=========================計算閾值========================
f=count（st+1:nd+1）;
p=st;q=nd-st;??????????????????????%p和q分別是灰度的起始和結(jié)束值
u=0;
for?i=1:q
????u=u+f（i）*（p+i-1）;
????ua（i）=u;
end
for?i=1:q
????w1（i）=sum（f（1:i））;
end
d=（u*w1-ua）.^2./（w1.*（1-w1））;???????%求出不同i值時類間方差
[ytp]=max（d）;??????????????????????%求出最大方差對應(yīng)的灰度級
th=tp+p????????????????????????????%threshold
if?th>200
????th=100;
end;
%?if?th<40
%?????th=45;
%?end;
W=zeros（ROWCOL）;???????????????????%新建一幅ROW×COL大小的圖

X=double（pic2）;

%%%從第2行第2列開始到倒數(shù)第2行倒數(shù)第2列

for?i=2:ROW-1

????for?j=2:COL-1
????????%%%%sobel算子
????????%%%%Dx是對X的偏導(dǎo)數(shù)，檢測水平邊緣，梯度方向是垂直的.
????????Dx=X（i+1j-1）+2*X（i+1j）+X（i+1j+1）-X（i-1j-1）-2*X（i-1j）-X（i-1j+1）;
????????%%%%Dy是對Y的偏導(dǎo)數(shù)，檢測垂直邊緣的梯度方向是水平的.
????????Dy=X（i-1j+1）+2*X（ij+1）+X（i+1j+1）-X（i-1j-1）-2*X（ij-1）-X（i+1j-1）;
????????%%%梯度的幅值??????
????????Mag=sqrt（Dx*Dx+Dy*Dy）;
????????%%%梯度的方向角Dy有取0的可能.
????????%%%a（ij）=round（180*（atan（Dx/Dy）））;
????????%%%梯度取整賦值給輸出圖像Wanted=W;
????????W（ij）=round（Mag）;
????end
end

for?i=1:ROW
????for?j=1:COL
????????if?W（ij）>th
???????????W（ij）=255;
????????else
???????????W（ij）=0;
????????end
????end
end
subplot（232）;
imshow（W）;
title（‘Sobel算子邊緣提取‘）;

for?i=2:ROW-1

????for?j=2:COL-1
????????%%%%prewitt算子
????????Dx=X（i+1j-1）+X（i+1j）+X（i+1j+1）-X（i-1j-1）-X（i-1j）-X（i-1j+1）;
????????Dx=Dx*4/3;
????????Dy=X（i-1j+1）+X（ij+1）+X（i+1j+1）-X（i-1j-1）-X（ij-1）-X（i+1j-1）;??
????????Dy=Dy*4/3;
????????Mag=sqrt（Dx*Dx+Dy*Dy）;
????????%%%梯度的方向角Dy有取0的可能.
????????%%%a（ij）=round（180*（atan（Dx/Dy）））;
????????%%%梯度取整賦值給輸出圖像Wanted=W;
????????W（ij）=round（Mag）;
????end
end

%%%閾值人為設(shè)定隨不同圖像不同目的改變.
%%%大于閾值認為是邊緣設(shè)為最亮;小于則為背景最黑

for?i=1:ROW
????for?j=1:COL
????????if?W（ij）>th
???????????W（ij）=255;
????????else
???????????W（ij）=0;
????????end
????end
end
subplot（233）;
imshow（W）;
title（‘Prewitt算子邊緣提取‘）;

for?i=1:ROW-1

????for?j=1:COL-1
????????%%%%prewitt算子
????????Dx=X（ij