利用迭代法求閾值並實現二值化

時間 2021-01-13

場景設定

將閾值計算的迭代法，設計爲函數 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*hi/w(T0)$
$u2=\sum_{i=T0}^{255}i*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;