利用迭代法求阈值并实现二值化

场景设定

将阈值计算的迭代法，设计为函数 level = thresh_x( f ); 并调用函数测试：读入lena.ppm，lena1.ppm, lena2.ppm, ocr.ppm等测试。

分析

这四个测试点不要想的太简单，图像ocr.ppm和图像lena.ppm这两个是最基本的测试点，只要你迭代法正确编写，就可以得出答案，但是lena1.ppm图像是整体偏暗，而lena2.ppm整体偏亮，这就意味着你迭代法的初始值需要一个算法来进行设定。这里仅仅是大致说说，听不懂没关系，只要知道一句话就好了，那就是没那么简单，我会一步一步来讲。(*^▽*)

迭代法思路

思路一（正解）

① 先图中灰度的中值作为初始阈值T0
② 利用阈值T0将图分割为两个部分（小于灰度值T0和大于灰度值T0），计算着两部分的灰度均值u1和u2
$$u1=\sum _{i=0}^{T0}i\ast hi/w(T0)$$
$$u2=\sum _{i=T0}^{255}i\ast hi/(1-w(T0))$$
其中$h(i)$是归一直方图，而$w(T0)$是归一累积直方图
③ 计算新的阈值 $T1=1/(u1+u2)$
④ 重复②和③步骤，直至T0和T1的差值小于某个给定值
⑤ 然后这个T0就是阈值了

代码实现

%灰度直方图
h=imhist(f);

%累积直方图
vmax = 256 ; H(1:vmax)=0; H(1)=h(1);
for i=2:vmax
    H(i) = H(i-1) + h(i); 
end;

% 归一化直方图
h1 = h / (sx * sy);

%归一累积直方图
H1 = H/(sx*sy);

while 1
    FZ1=0;FZ2=0; %分子
    for i=1:256
        if i<T0
           FZ1=FZ1+(i-1)*h1(i);
        else
           FZ2=FZ2+(i-1)*h1(i);
        end;
    end;
    u1=FZ1/H1(T0);
    u2=FZ2/(1-H1(T0));
    T1=(u1 + u2)/2;
    if abs(T0-T1)<1
        %求出灰度值比例
        level = double(T0);
        level = level / 255;
        break;
    else
        T0=abs(uint8(T1));
    end;
end;

这个思路有点绕，毕竟用到了归一直方图和归一累积直方图的概念，需要好好去理解下，下面这种思路就比较的简单，完全初中知识，应该好理解一点~~

思路二

首先我们知道要算一个整体图像的像素平均值，就等于用该点像素值乘以值为该像素值的像素个数，然后相加起来，最后除以总共的像素个数，这就算出了整体图像的像素平均值，这里就是利用这个思路。

代码实现

while 1
    u1=0; u2=0;
    FM1=0;FM2=0;
    for i=1:255
        u1_ave=0; u2_ave=0;
        if i<T0
            u1=u1+h(i)*(i-1);
            FM1=FM1+h(i);
            u1_ave=u1/FM1;
        else
            u2=u2+h(i)*(i-1);
            FM2=FM2+h(i);
            u2_ave=u2/FM2;
        end;
        T1=(u1_ave+ u2_ave)/2;
    end;
    if T0-T1<1
        break;
    else
        T0=T1;
    end;
end;

我以下使用的都是思路一的方式来求的，原因是思路二好像是超过int型能表示的最大数，导致弄出来的图像效果没有思路一的好些~~

如果我们理解了迭代法的思路后，我们就可以编写出该函数 level = thresh_x( f )，来试试了。

测试

thresh_x.m

function level = thresh_x(f)
[sx,sy]=size(f);

T0=30; %直接设置初始值

%灰度直方图
h=imhist(f);

%累积直方图
vmax = 256 ; H(1:vmax)=0; H(1)=h(1);
for i=2:vmax
    H(i) = H(i-1) + h(i); 
end;

% 归一化直方图
h1 = h / (sx * sy);

%归一累积直方图
H1 = H/(sx*sy);

while 1
    FZ1=0;FZ2=0; %分子
    for i=1:256
        if i<T0
           FZ1=FZ1+(i-1)*h1(i);
        else
           FZ2=FZ2+(i-1)*h1(i);
        end;
    end;
    u1=FZ1/H1(T0);
    u2=FZ2/(1-H1(T0));
    T1=(u1 + u2)/2;
    if abs(T0-T1)<1
        %求出灰度值比例
        level = double(T0);
        level = level / 255;
        break;
    else
        T0=abs(uint8(T1));
    end;
end;

函数编写完了，我们就编写一个thresh_x_test.m来测试我们的函数把

close all; clear all;

f= imread('../images/ocr.ppm');
level = thresh_x(f);
g = im2bw(f,level);
figure;imshow(g);title(['res',num2str(level)]);

效果

但是在弄lena1.ppm和lena2.ppm的时候就一直在报错，以下就以lena2.ppm来说明问题

我们通过设置断点可以发现，在第一次计算u1的时候，它的值为NAN，为啥会这样呢，通过分析可以发现我们在求u1和u2的时候是使用归一/归一累积，而第一次的T0是我们自己指定的，我指定为30，但是呢，由于lena2.ppm整体图像偏亮，导致在前100的像素占比为0，那么在算$u1=FZ1/H1(T0)$时，H1(T0)为0，所以u1就会变成NAN了 （如果我们把T0指定为150左右则可以正常弄lena2.ppm，但是在弄lena1.ppm时就需要把T0指定为35左右）。所以我们就需要去写一个算法，让它自动判断大致的范围。

我们最自然的想法是直接找出一个图像中的最大像素和最小像素，然后求一次平均就可以了，代码如下所示：

% 初始值设定
ZMAX=max(max(f));
ZMIN=min(min(f));
T0=(ZMAX+ZMIN)/2;

这样确实能处理大多数图像，但是存在极端情况，就会让这种方法失效。举个栗子来说，假如一张图像整体偏亮，大部分像素都分布在200以上，但是呢，就是有一个像素，它的像素值就是0，根据这个算法，它是找整张图的最大像素值和最小像素值来求平均，那么可能求出的这个平均的像素值可能是128，那么此时，128对应的归一累积就会很小，在求u1时，就会无穷大出现NAN的情况。所以这个算法并不合适。
为了避免这种情况，我的思路是求归一累积为0.05左右的和归一累积为0.95左右对应的灰度值，用这两个灰度值求平均，代码如下：

%求初始值
H2=roundn(H1,-2);
for i=1:256
    if 0.03<H2(i) && H2(i)<0.06
        ZMIN=i;
    end;
    if 0.95<H2(i) && H2(i)<0.98
        ZMAX=i;
    end;
end;
T0=(ZMIN+ZMAX-2)/2;
T0=fix(T0); %取整

使用这个方法来大致获取一个合适的初值就可以避免极端情况了。

整体代码

function level = thresh_x(f)

[sx,sy]=size(f);
%灰度直方图
h=imhist(f);

%累积直方图
vmax = 256 ; H(1:vmax)=0; H(1)=h(1);
for i=2:vmax
    H(i) = H(i-1) + h(i); 
end;

% 归一化直方图
h1 = h / (sx * sy);

%归一累积直方图
H1 = H/(sx*sy);

%求初始值
H2=roundn(H1,-2);
for i=1:256
    if 0.03<H2(i) && H2(i)<0.06
        ZMIN=i;
    end;
    if 0.95<H2(i) && H2(i)<0.98
        ZMAX=i;
    end;
end;
T0=(ZMIN+ZMAX-2)/2;
T0=fix(T0); %取整

%迭代法
while 1
    FZ1=0;FZ2=0; %分子
    for i=1:256
        if i<T0
           FZ1=FZ1+(i-1)*h1(i);
        else
           FZ2=FZ2+(i-1)*h1(i);
        end;
    end;
    u1=FZ1/H1(T0);
    u2=FZ2/(1-H1(T0));
    T1=(u1 + u2)/2;
    if abs(T0-T1)<1
        %求出灰度值比例
        level = double(T0);
        level = level / 255;
        break;
    else
        T0=abs(uint8(T1));
    end;
end;

效果