图像二值化的阈值求法

图像二值化的操作，就是寻找一个合适的阈值T，把图像中所有像素值小于等于该阈值的像素点设置为一个像素值，同时把图像中所有像素值大于该阈值的像素点设置为另一个像素值。也就是说，二值化之后整幅图像中只有两种灰度值。对于一帧8位的图像，二值化操作通常是把其像素值设置为0或255，从直观上看就是把具有灰度变化的图像设置为只有黑与白的黑白图像，0值对应黑色，255值对应白色。由此可知，二值化效果的好坏取决于阈值T的选取是否合适，阈值T就相当于黑白色的分界线。

图像二值化的操作在数字图像处理中具有重要作用，该操作把图像变成黑白之后，可以凸显图像的重要轮廓与结构，大大减少运算量。比如在细胞的荧光成像中，常常需要把细胞区域（也即较亮区域）识别出来，很重要的一步就是把荧光成像的图像转换为二值图像，细胞区域为白色，背景区域为黑色，这样才便于检测细胞轮廓，从而识别细胞。

最常用的求阈值算法有大津法和迭代法。本文针对荧光图像的二值化操作，详细讲解这两种算法的原理，同时对迭代法进行改进，使其可以求出更加合适的阈值。

1. 大津法，简称为OTSU算法，属于自适应的阈值确定方法，其核心思路是寻找一个阈值T，把图像的所有像素点分成两类，一类的像素值均小于等于T（背景区域），另一类的像素值均大于T（前景区域），当这两类的类间方差取得最大值时，则认为该T值为最合适的阈值。假设背景区域所有像素点数为n0，其占图像总像素数的比例为p0，平均像素值为m0；前景区域所有像素点数为n1，其占图像总像素数的比例为p1，平均像素值为m1；图像的所有像素点数为n、平均像素值为m。那么有以下关系式：

定义类间方差g为：

那么由以上3式可以得到下式，这就是类间方差的计算公式。

循环遍历图像可能取到的所有像素值作为阈值T，使g取得最大的T值则认为是最合适的阈值。比如8位图像所有可能的取值为0~255，那么分别计算出T取值0~255时的g值，当g值取得最大值，则认为找到了最佳阈值。

使用大津法确定阈值，对一张荧光细胞图像进行二值化，得到的结果如下图所示。可以看到细胞区域基本被标识为白色区域，与背景区域区别开来。

原图

大津法求阈值的二值化结果图

2. 迭代法

迭代法也属于自适应的阈值确定方法，其核心思路是通过多轮迭代运算来一步一步逼近最佳阈值。设置一个初始阈值T0，同样把图像的所有像素点分成两类，一类的像素值均小于等于T0（背景区域），另一类的像素值均大于T0（前景区域），然后分别计算背景区域与前景区域的像素平均值m0、m1，并取m0与m1之和的一半作为下一轮迭代的阈值。如果当前轮迭代得到的阈值与上一轮迭代的阈值相差很小，则结束迭代。迭代流程如下图所示：

使用迭代法求阈值，并对同一张荧光细胞图像进行二值化，得到的结果如下图所示。可以看到细胞区域也基本被识别出来。

迭代法求阈值的二值化结果图

3. 迭代法的改进。当荧光图像的背景亮度与细胞区域亮度的区别不是那么明显，也即背景区域比较亮的时候，使用大津法和迭代法确定阈值，会得到比较差的二值化结果，细胞区域不能很好地与背景区域区别开来。如下图所示，可以看到很大片的背景区域二值化之后也变成了白色。

原图

大津法求阈值的二值化结果图

迭代法求阈值的二值化结果图

由上图可知，由大津法与迭代法取得的阈值，与较亮的背景区域的像素值比较接近，导致较亮的背景区域也被误划分为细胞区域（白色区域）。下面我们对迭代法进行改进，使其求得的阈值更接近细胞区域的像素值，从而与较亮背景区域的像素值区分开来。观察迭代法的阈值更新公式，可以知道该算法每次取m0与m1之和的一半来更新阈值：

为了使阈值更接近m1，我们修改上述公式，增大m1的比重，同时减小m0的比重：

matlab的graythresh函数已经实现了大津法求阈值，所以在此我们不再重复实现，下面仅贴出改进之后的迭代法matlab代码：

% 最佳阈值法获取图像二值化的阈值
% T - 得到的最佳阈值
% I - 输入图像
% allow - 迭代过程中，当前迭代得到的阈值与下一轮迭代得到的阈值之差的绝对值小于allow，则认为找到了最佳阈值


function [T] = get_threld(I, allow)


[x,y] = size(I);                % 求出图象行列数
b = double(I);       %把图像由整型转换为浮点型           
                        
T_current = mean(b(:));  %取图像的均值作为初始阈值T0


T_next = 0;                               %TT赋初值
S0 = 0.0;               %为计算灰度大于阈值的元素的灰度总值、个数赋值
n0 = 0.0;                       
S1 = 0.0;                  %为计算灰度小于阈值的元素的灰度总值、个数赋值
n1 = 0.0;                      
d = abs(T_current-T_next);
count = 0;
% 记录几次循环
while(d >= allow)                 % 迭代最佳阈值分割算法
    
    count = count+1;       %迭代次数加1
    
    for i = 1:x      %遍历每一行
        for j = 1:y  %遍历每一列
            
            if b(i,j) >= T_current   %计算图像中所有像素值大于阈值的点的像素值之和
                S0 = S0 + b(i,j);     
                n0 = n0 + 1;
            else             %计算图像中所有像素值小于阈值的点的像素值之和
                S1 = S1 + b(i,j);
                n1 = n1 + 1;
            end


        end
    end 
    
    T0 = S0/n0;   %计算图像中所有像素值大于阈值的点的像素值平均值
    T1 = S1/n1;   %计算图像中所有像素值小于阈值的点的像素值平均值
    
    %T_next = (T0+T1)/2;     %原迭代算法
    T_next = T1 + (T0-T1)*0.9;  %改进的迭代算法
    
    d = abs(T_current-T_next);    %当前迭代得到的阈值与下一轮迭代得到的阈值之差的绝对值
    T_current = T_next;         %把下一轮迭代的阈值赋值给当前轮迭代的阈值
end


mmax = max(max(I));
T = T_current/double(mmax);       %把求得的阈值处以图像的最大像素值进行归一化


end

运行以上matlab代码，得到的二值化结果如下图所示。可以看到改进之后的迭代算法求取的阈值，可以很好地区分细胞区域与背景区域。因此我们的小改进还是具有不错的效果的。

改进的迭代法求阈值的二值化结果图