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

场景设定

将阈值计算的迭代法,设计为函数 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
u 1 = ∑ i = 0 T 0 i ∗ h i / w ( T 0 ) u1=\sum_{i=0}^{T0}i*hi/w(T0) u1=i=0T0ihi/w(T0)
u 2 = ∑ i = T 0 255 i ∗ h i / ( 1 − w ( T 0 ) ) u2=\sum_{i=T0}^{255}i*hi/(1-w(T0)) u2=i=T0255ihi/(1w(T0))
其中 h ( i ) h(i) h(i)是归一直方图,而 w ( T 0 ) w(T0) w(T0)是归一累积直方图
③ 计算新的阈值 T 1 = 1 / ( u 1 + u 2 ) T1=1/(u1+u2) 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,那么在算 u 1 = F Z 1 / H 1 ( T 0 ) u1=FZ1/H1(T0) 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;
效果

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

    原文作者:duganlx
    原文地址: https://blog.csdn.net/qq_40626497/article/details/89702723
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系博主进行删除。
点赞