[數字圖像處理]灰階變換——直方圖處理

直方圖均衡

        一副圖像的直方圖，表示了其灰階分布的特性。對于數字圖像來說，假設灰階值k出現了

次，那麼其機率密度函數如下所示。

這個式子，表示了像素的灰階值為k機率。其中，M與N為圖像的尺寸。

       對于一幅動态範圍較窄的圖像，其歸一化灰階直方圖如下所示。

       對于此類圖像，使用灰階拉伸，也能使其的動态範圍得到改善，增強圖像對比度。但是，灰階拉伸從本質上來講，是很暧昧的，灰階拉伸沒有什麼明确的目的，隻是依靠某個函數的做了灰階變換，使得圖像的在灰階直方圖的分布範圍擴大。灰階拉伸對于結果沒有嚴格的要求，是以根據變換函數的不同，灰階拉伸有無數種結果。

       與灰階拉伸不同的是，直方圖均衡的目的就比較明确，使用某個特殊的函數，使原圖像的灰階分布平均化。為了将其平均化，這裡需要用到累積分布函數（Probability Density Function）的概念。

所謂的平均化，直方圖均衡的目的是将一幅圖的累積分布的曲線（下圖左），變為下圖右的分布形式。

                                                 （左）原圖的累積分布曲線                        （右）理想化直方圖均衡後的累積分布曲線

       從右圖中可以看出，為了保證直方圖均衡後的累積分布曲線是一條傾角為45°的直線，那麼，變換後的函數的機率密度函數應與輸入值無關，且保持一個常量。輸入r與輸出s有如下關系，

将變換關系

帶入上式子，可得

将其帶入輸入r與輸出s關系，可得到

到此，已經可以看出，輸出的機率密度函數，與s的取值無關。很明顯，這個分布是絕對平均的。

       但是，事實上，上面的式子均是在連續的情況下的。在數字圖像的領域，我們必須将其離散化，才能使用。離散的情況下，其累計分布函數如下所示。

至此，直方圖均衡就完成了。

       為了驗證結果，其直方圖均衡後的圖像與灰階直方圖如下所示。

結果就是這樣了。但是，呃！！！！問題來了，這是直方圖均衡麼？這跟灰階拉升有何差別，而且其均衡後的圖像的直方圖也不是一個很平均的值，和我們的構想有很大差別。這個結果真的沒問麼？為了回答這個問題，将其原圖像的累積分布函數與直方均衡後的累積分布函數畫出來，如下。

                                                      （左）原圖的累積分布曲線                                             （右）直方圖均衡後的累積分布曲線

從累積分布曲線中，我們可以看到，直方圖均衡後的曲線，确實跟我們設想的理想曲線很接近！由于我們證明的時候，是在連續的情況下去思考的，是以，将其離散之後，多少是有不理想的情況發生的。均衡後的累積分布曲線确實是一個不錯的結果。這樣看來，我們才能看清楚直方圖均衡的本質。若不是這樣，還真不知道直方均衡與灰階拉伸的差別。

為此，直方圖均衡成功，Matlab代碼如下。

close all;
clear all;

%% -------------Histogram Equalization-----------------
close all;
clear all;

f = imread('washed_out_pollen_image.tif');
f = mat2gray(f,[0 255]);

[M,N] = size(f);
g = zeros(M,N);

r = imhist(f)/(M*N);  %(500*500 is size of image)
s = zeros(256,1);

for k = 1:1:256
    for j = 1:1:k
        s(k,1) = s(k,1) + r(j);
    end
end

for x = 1:1:M        %(500*500 is size of image)
    for y = 1:1:N
        g(x,y) = s(uint8(f(x,y)*255)+1,1);
    end
end 


figure();
subplot(2,2,1);
imshow(f);
xlabel('a).Original Image');

subplot(2,2,2);
h = imhist(f)/(M*N);
bar(0:1/255:1,h);
axis([0 1 0 0.1]),grid;
%axis square;
xlabel('b).The Histogram of a');
ylabel('Number of pixels');

subplot(2,2,3);
imshow(g);
xlabel('c).Histogram Equalization');

subplot(2,2,4);
h = imhist(g)/(M*N);
bar(0:1/255:1,h);
axis([0 1 0 0.1]),grid;
%axis square;
xlabel('d).The Histogram of c');
ylabel('Number of pixels');

figure();
x=0:1/255:1; 
plot(x,s);
axis([0,1,0,1]),grid;
axis square;
xlabel('intensity level of r');
ylabel('P_{r}(r)');

r = imhist(g)/(M*N);
s = zeros(256,1);
for k = 1:1:256
    for j = 1:1:k
        s(k,1) = s(k,1) + r(j);
    end
end
figure();
x=0:1/255:1; 
plot(x,s);
axis([0,1,0,1]),grid;
axis square;
xlabel('intensity level of s');
ylabel('P_{s}(s)');
           

（作者注：其實，直方圖均衡有直接的函數，這裡隻是想按照自己了解的方法去做個均衡出來看看，是以進階一點的庫函數基本沒有，都用循環來的）

直方圖比對

       前面，我們說了直方圖的均衡，直方均衡可以得到一個灰階直方圖分布平均的圖像。但是有一些問題，直方圖均衡的結果是唯一的。這是其優點，也是其缺點。優點在于，無需多餘的參數，就可以實作對比度的增強。缺點在于，沒有參數，對于結果無法調整。舉個栗子（輸入法君賣萌了╭(╯3╰)╮），一幅很黑的圖像，使用直方圖均衡的話，這幅圖像的直方均衡結果，将使得圖像偏白。使用 《Digital Image Processing》 Rafael C. Gonzalez / Richard E. Woods 的例子來說，看下圖。

       如上圖所示，直方均衡的結果不是很理想。由直方均衡的結果看來，我們能看到，其實對比原圖，原圖一些細節已經顯露出來了，隻要稍加修改，就可以得到很理想的增強結果。為此，直覺告訴我們（不要吐槽這句話，有時候直覺就是那麼簡單。\(≧▽≦)/），隻需要将原圖的灰階直方圖稍微左移一些即可。但是，直方均衡卻是不可調整的，為此，就有了直方圖比對。

       對于原圖像，其機率分布函數是

其輸出也就是直方均衡的結果。為了實作直方圖比對，我們還需要一個函數，也就是我們所期待的機率分布函數。這個機率分布函數的機率分布函數如下所示。

其中，z代表了期望的分布。我們假設如下關系成立。

也就是，我們所期待的直方圖分布的機率分布函數與原圖像的直方圖機率分布函數相等。那麼，我們可以利用G()的反函數，我們也就可以得到我們所期待的直方圖。

       上面的描述可能太過拖沓，我們采用簡潔一點的描述。

       一. 将原圖進行直方圖均衡。

       二. 假設我們期待的分布，直方均衡的結果與原圖直方圖均衡的結果一緻。利用反函數，還原成我們期待的分布。

       所實作的直方圖均衡結果如下所示。

根據直方圖比對的結果，我們可以看出，所得出的結果和我們所期待的分布很接近。其結果也和我們的初衷一樣，使得原圖的直方圖稍微左移了一些。所得到的結果比起直方均衡要好得多。

其Matlab代碼如下。

%% -------------Histogram Matching-----------------
close all;
clear all;

f = imread('mars_moon_phobos.tif');
f = mat2gray(f,[0 255]);

[M,N] = size(f);
g = zeros(M,N);

r = imhist(f)/(M*N); 
s = zeros(256,1);

for k = 1:1:256
    for j = 1:1:k
        s(k,1) = s(k,1) + r(j);
    end
end
 
for x = 1:1:M        
    for y = 1:1:N
        g(x,y) = s(uint8(f(x,y)*255)+1,1);
    end
end 

figure();
subplot(2,2,1);
imshow(f);
xlabel('a).Original Image');

subplot(2,2,2);
h = imhist(f)/(M*N);
bar(0:1/255:1,h);
axis([0 1 0 0.1]),grid;
xlabel('b).The Histogram of a');

subplot(2,2,3);
imshow(g);
xlabel('c).Histogram Equalization');

subplot(2,2,4);
h = imhist(g)/(M*N);
bar(0:1/255:1,h);
axis([0 1 0 0.1]),grid;
xlabel('d).The Histogram of c');

%---------------pz-----------------------
r = r';
p = [0:r(1,1)/15:r(1,1) 17*r(1,2:241)];
%p = [0:r(1,1)/15:r(1,1) r(1,2:241)];
%p = [0:55:255 , 255:-25:24 , 24:-0.14:0 , 0:1:15 , 15:-0.285:0];
p = (p/sum(p));
p = p';
r = r';

G = zeros(256,1);
for k = 1:1:256
    for j = 1:1:k
        G(k,1) = G(k,1) + p(j);
    end
end


G_Negatives = zeros(256,1);
Flag = 0;
for k = 1:1:256
    for j = 1:1:256
            if(uint8(G(j,1)*255) == k) 
                 G_Negatives(k,1) = (j-1)/255;
                 Flag = 1;
            end
    end
    
    if(Flag == 0)
        if(k == 1) G_Negatives(k,1) = 0;
        else G_Negatives(k,1) = G_Negatives(k-1,1); end   
    end
    Flag = 0;

end

z = zeros(M,N);
for x = 1:1:M       
    for y = 1:1:N
        z(x,y) =  G_Negatives(uint8(g(x,y)*255)+1,1);
    end
end 

figure();
subplot(2,2,1);
imshow(z);
xlabel('a).Histogram Matching');

subplot(2,2,2);
h = imhist(z)/(M*N);
bar(0:1/255:1,h);
axis([0 1 0 .1]),grid;
xlabel('b).The Histogram of a');

subplot(2,2,3);
x=0:1/255:1; 
plot(x,p);
axis([0,1,0,.1]),grid;
xlabel('c).Specifieds Histogram (Normalized)');
 
subplot(2,2,4);
x=0:1/255:1; 
plot(x,G,x,G_Negatives);
axis([0,1,0,1]),grid;
axis square;
xlabel('Input intensity level');
ylabel('Onput intensity level');
           

部落格位址：http://blog.csdn.net/thnh169/