基于阈值的灰度图像提取法

　　对于简单的灰度图像，如果目标与背景的灰度存在一定的差异，那么可以用阈值来提取目标。关键是确定阈值，常用方法有：

模态法、P参数法、可变阈值法、大津法和迭代逼近法等。

模态法：

　　取直方图的波谷作为阈值。适用于目标与背景灰度差异较大，目标与背景的直方图各有一个波峰的灰度图像。

　　如果直方图凹凸较激烈，寻找波谷存在困难，可以用领域点平均化处理使直方图相对平滑。

P参数法：

　　当目标占图像的面积比例已知，且目标最小灰度值大于背景最大灰度值时，可采用P参数法。具体做法是，如果目标灰度接近白色，则在

　　直方图右侧开始累计像素，直到累计像素占总像素数比例为P时，将当前的灰度作为阈值（否则从左侧开始累计）。

　　例如下图，月亮所占比列为0.45，它最低灰度与背景灰度稍有重合，阈值确定为75，最终基本提取出目标。

　　　　　　　　

可变阈值法：

　　当图像背景灰度多变时，可以为图像不同部位设置不同的阈值，这样可以很好地去除干扰。

　　例如提取下图的蒲公英，图像上部光照多，主要为蒲公英的羽毛，下部光照少，且主要为茎，对上中下采取不同阈值的提取效果会好很多。

　　

　　　　　　　　　　

大津法：

　　大津法的应用最广泛，有很好的自适应阈值。适用于批量图像提取或者动态灰度图像提取（然而事实证明，当大津遇到反光时，干扰也会很严重）。

　　大津法的原理如下：

　　设T为分割阈值，由此划分出区域记为R₁和R₂。

　　R₁占图像总像素的比例为θ₁，R₂占图像总像素的比例为θ₂。

　　图像的平均灰度值为μ₀，R₁的平均灰度为μ₁，R₂的平均灰度为μ₂。

　　显然有：

　　　　μ₀=θ₁*μ₁+θ₂*μ₂

　　如果阈值选得恰当，则R₁与R₂的灰度差应该较大，那么μ₁与μ₀、μ₂与μ₀的差的绝对值也就较大。

　　定义类间方差，来反应灰度差大小：

　　　　σ²=θ₁（μ₁-μ₀）²+θ₂（μ₂-μ₀）²

　　对于每个选定的阈值T，都有σ对应，而最大σ的阈值T则是最佳阈值。

　　另外，尽管有化简公式来减少迭代次数的方法（https://www.cnblogs.com/kensporger/p/11270452.html），大津法的过程计算量还是很大的。

迭代逼近法：

　　此法是大津法的简化版，具体步骤如下：

　　取最大灰度与最小灰度平均值（也可以是图像平均灰度）为初始阈值T，由此得两个区域为R₁和R₂

　　计算R₁和R₂的灰度均值μ₁、μ₂

　　选择新阈值T=（μ₁+μ₂）/2，重复以上步骤,直到T不再变化或变化在一定范围内。

　　采用该法的效果如图所示，大概只迭代了五次左右：

　　

　　但是，该法确定的阈值很粗糙，往往与最佳阈值会有较大偏差，比如提取之前的月亮，表现得就并不是很好了：

　　

以上大部分matlab仿真测试代码罗列如下：

%p参数法确定阈值
function thres  =   pthres(file,p)

    imga = imread(file);
    
    %灰度化
    imga = rgb2gray(imga);
   
    %获取大小
    [sizex,sizey]=size(imga);
    
    %直方图
    histogram = zeros(256,1);
    for i=1:sizex
        for j=1:sizey
            histogram(imga(i,j)+1)=histogram(imga(i,j)+1)+1;   
        end
    end
    
    %累计像素
    pixels_cnt=0;
    thres=256;
    while(pixels_cnt<p*sizex*sizey)
        pixels_cnt=pixels_cnt+histogram(thres);
        thres=thres-1;
    end
    
    %根据阈值二值化
    mybw=imbinarize(imga,(thres-1)/255);
    
    subplot(1,2,1);
    imshow(imga);
    title('原图');
    subplot(1,2,2);
    imshow(mybw);
    title('p参数法二值化');
end

%多阈值分割实例
imga=imread('multi_thres.jpg');
imga=rgb2gray(imga);

single  = imbinarize(a,85/255);
multi(250:300,:) = imbinarize(a(250:300,:),85/255);
multi(150:250,:) = imbinarize(a(150:250,:),150/255);
multi(1:150,:) = imbinarize(a(1:150,:),170/255);

subplot(1,3,1);
imshow(imga);
title('原图');
subplot(1,3,2);
imshow(single);
title('单阈值处理');
subplot(1,3,3);
imshow(multi);
title('多阈值处理');

%迭代逼近法
function thres=iterator(file)
    imga = imread(file);
    
    %灰度化
    imga = rgb2gray(imga);
   
    %获取大小
    [sizex,sizey]=size(imga);
    
    %直方图与总灰度统计
    gray_sum=0;
    histogram = zeros(256,1);
    for i=1:sizex
        for j=1:sizey
            histogram(imga(i,j)+1)=histogram(imga(i,j)+1)+1; 
            gray_sum= gray_sum+double(imga(i,j)+1);
        end
    end
    %初始阈值
    thres=gray_sum/(sizex*sizey);
    old_thres=0;
    %迭代逼近
    while(abs(thres-old_thres)>1)
       disp(thres);
        u1=0;u2=0;
        cnt1=0;cnt2=0;
        for i=1:thres
            u1=u1+histogram(i)*i;           %统计R1区域平均灰度
            cnt1=cnt1+histogram(i);
        end
        for i=256:-1:thres
            u2=u2+histogram(i)*i;           %统计R2区域平均灰度
            cnt2=cnt2+histogram(i);
        end   

        u1=u1/cnt1;
        u2=u2/cnt2;

        old_thres=thres;
        thres=(u1+u2)/2;                    %新阈值
    end
    
    mybw=imbinarize(imga,(thres-1)/256);

    subplot(1,2,1);
    imshow(imga);
    title('原图');
    subplot(1,2,2);
    imshow(mybw);
    title('迭代逼近法二值化');    
    
end