前几天看论文用到阈值选取算法，现用matlab实现OSTU算法如下：

  Otsu最大类间方差法原理

    利用阈值将原图像分成前景，背景两个图象。

    当取最佳阈值时，背景应该与前景差别最大，即方差最大。otsu算法找的就是这个最大方差下的阈值。

 最大类间方差法（otsu）的公式推导：

   记t为前景与背景的分割阈值，前景点数占图像比例为w0， 平均灰度为u0；背景点数占图像比例为w1，平均灰度为u1。

   则图像的总平均灰度为：u=w0*u0+w1*u1。

   前景和背景图象的方差：g=w0*(u0-u)*(u0-u)+w1*(u1-u)*(u1-u)=w0*w1*(u0-u1)*(u0-u1),此公式为方差公式。

   循环求取最大方差即可。

   MABLAB代码及详细注释：

 

 

function  ostu

img = imread('Lena.jpg');

I_gray=rgb2gray(img);

 

figure,imshow(I_gray);

I_double=double(I_gray);    %转化为双精度，因为大多数的函数和操作都是基于double的

%figure,imshow(I_double);

[wid,len]=size(I_gray);     %wid为行数，len为列数

colorlevel=256;     %灰度级

hist=zeros(colorlevel,1);   %直方图，256×1的0矩阵

 

%threshold=128; %初始阈值

%计算直方图，统计灰度值的个数

for i=1:wid

    for j=1:len

        m=I_gray(i,j)+1;    %因为灰度为0-255所以+1

        hist(m)=hist(m)+1;

    end

end

 

%直方图归一化

hist=hist/(wid*len);   

 

%miuT为总的平均灰度，hist[m]代表像素值为m的点个数

miuT=0;

for m=1:colorlevel

    miuT=miuT+(m-1)*hist(m);

end

xigmaB2=0;  %用于保存每次计算的方差，与下次计算的方差比较大小

for mindex=1:colorlevel

    threshold=mindex-1;

    omega1=0;    %前景点所占比例

    omega2=0;    %背景点所占比例

    for m=1:threshold-1

        omega1=omega1+hist(m);  %计算前景比例

    end

    omega2=1-omega1;            %计算背景比例

    miu1=0;      %前景平均灰度比例

    miu2=0;      %背景平均灰度比例

    %计算前景与背景平均灰度

    for m=1:colorlevel

        if m

        miu1=miu1+(m-1)*hist(m);   %前景平均灰度

        else

        miu2=miu2+(m-1)*hist(m);   %背景平均灰度

        end

    end

   

 %   miu1=miu1/omega1;

%    miu2=miu2/omega2;

    %计算前景与背景图像的方差

    xigmaB21=omega1*(miu1-miuT)^2+omega2*(miu2-miuT)^2;

    

    xigma(mindex)=xigmaB21;   %保留每次计算的方差

    

     %每次计算方差后，与上一次计算的方差比较。保留最大方差为 finalT

    if xigmaB21>xigmaB2

        finalT=threshold;

        xigmaB2=xigmaB21;

    end

end

 

%比较方法两种阈值的不同

fT=finalT/255;         %阈值归一化

T=graythresh(I_gray);  %matlab函数求阈值

 

for i=1:wid

    for j=1:len

        if I_double(i,j)>finalT

            bin(i,j)=255;

        else

            bin(i,j)=0;

        end

    end

end

 

figure,imshow(bin);

figure,plot(1:colorlevel,xigma);

 

end

 

运行结果：

http://s6/mw690/a98e39a2gcc18aa1a0e15&690