數字影像直方圖匹配或規定化Histogram Matching (Specification)處理

? ? 前往老猿Python博文目錄 https://blog.csdn.net/LaoYuanPython ?

一、引言

在《數字影像處理：直方圖均衡（Histogram Equalization）的原理及處理介紹》了直方圖均衡，直方圖均衡能自動地確定均衡需要的變換函式，該函式能產生有灰度值均勻直方圖的輸出影像，當需要對一副影像自動增強時，這是一種好方法，因為這種處理結果可以預知（predictable）、實作簡單，但是對某些應用，基于統一的直方圖的基礎增強（base enhancement）并不是最好的方法，特別地有時處理人員希望能指定結果影像的直方圖形狀（specify the shape of the histogram），這時就需要用到直方圖匹配（規定化），

本文是《數字影像處理》第三章直方圖處理中的直方圖匹配或規定化相關章節的內容，閱讀本文需要具備導數和微積分相關知識，如果需要復習或學習相關高數知識，請參考《人工智能數學基礎》，

二、知識概要：直方圖匹配或規定化

2.1、概述

直方圖匹配（Histogram Matching ），又稱為直方圖規定化（Histogram Specification），是指在處理影像時，由處理人員指定輸出影像的直方圖形狀，

老猿注：

• 指定直方圖形狀應該是通過指定直方圖的PDF函式來實作的；
• 本部分的原理沒有采用原文相同的方式介紹，因為原文比較拗口，乍一看有點難以理解，老猿按自己的理解來介紹，

2.2、連續灰度值的非數字影像直方圖匹配或規定化

2.2.1、相關原理

在《數字影像處理：直方圖均衡（Histogram Equalization）的原理及處理介紹》介紹了對于非數字影像的灰度值連續影像，對于灰度值為r的影像到直方圖均衡后灰度為s的影像有如下公式：

現假設針對該輸入影像，存在另一個灰度值集合z的影像，其影像內容與灰度值集合為r的影像是同一副影像，只是灰度值不同（這里的r、z、s在公式內是表示單個像素的灰度值，適用于影像的每個像素，同時在描述時r、z、s又表示了同一個影像所有像素的灰度集合的陣列），

類似式（3.3-10），z對應的灰度影像通過如下函式G（z）可以轉換為s對應的影像：

即可得：

這是因為針對同一影像有不同灰度分布的情況，都可以通過直方圖均衡化為灰度值均衡的影像，這就意味著兩副不同灰度值的同樣內容的影像都可以通過直方圖均衡后的結果影像關聯起來，并以此為中介求得兩副相同內容但灰度值不同的影像之間的灰度的映射關系，

由式(3.3-12) z = G^-1|T（r）|可知：z和r的映射只需要z和r對應的兩副影像的灰度值的概率密度函式經過定積分以及四則運算就可以求出，

這就是說，兩副內容相同灰度值不同的影像，可以很方便地找到兩副影像的灰度映射關系，這也就是為什么該灰度變換被稱為直方圖匹配或直方圖規定化的原因，它們的目的就是給定兩副相同內容不同灰度分布影像，通過這兩副影像本身的灰度概率密度函式就可以找到二者的映射關系，

2.2.2、連續灰度影像的直方圖規定化處理步驟

通過以上三個式子可知，使用下列步驟，可由一幅給定影像得到一幅其灰度級具有指定概率密度函式的影像：

1. 首先由輸入影像的概率密度函式PDF使用式(3.3-10)對輸入影像進行均衡得到像素值為s的輸出影像；
2. 使用式(3.3-11)結合輸出影像的PDF求得輸出影像到直方圖均衡像素值為s的影像的映射函式G（z）;
3. 對均衡后的影像像素值s進行G（z）的反變換G^-1（s），就可以得到輸出影像的像素值，

當所有像素都轉換后，就可以得到PDF等于指定函式的輸出影像，

2.2.3、連續灰度影像的直方圖規定化案例

例：假設采用連續灰度值的一幅影像的灰度 PDF 為pr（r）= 2r/(L-1)2，0≤r≤(L-1)，對于其他r值有pr（r）=0，尋找一個變換函式，使得產生的影像的灰度PDF是pz(z)=3z2/(L-1)3，0≤z≤(L-1)，而對于其他z值有pz(z)=0，

注：上面pr、pz的字母r和z都是下標，表示不同的概率密度函式，為了輸入方便就寫成了pr、pz，下同，

解：
首先，我們對區間[0，L-1]尋找直方圖均衡變換：

由定義可知，對于范圍[0，L-1]外的值，該變換為0，按直方圖均衡變換的處理可知，輸人灰度值的平方除以(L-1)（老猿注：原文和譯本這里都說是（L-1）2，經筆者仔細核實應該沒有平方）將產生一幅灰度為s并具有均勻PDF的影像，

為了得到pz(z)=3z2/(L-1)3指定直方圖的輸出影像，在[0，L-1]區間上，pz（z）對應的直方圖均衡變換公式為：

由定義可知在該區間外這個函式也為0，最后，我們要求G(z)=s. 但G(z)=z3/(L-1)2；因此z3/(L-1)2=s，可得：

這表明如果我們用(L-1)2乘以直方圖均衡過的每一個像素，取該乘積的 1/3 次冪，結果將是一幅我們所期望的影像，該影像在區間[0，L-1]內z的灰度的PDF為pz(z)=3z2/(L-1)3，

因為s=r2/(L-1)，我們可以直接由輸入影像的灰度r生成z：

這樣原影像中每一個像素值的平方與(L-1)相乘，然后再取該乘積的1/3次幕，將得到其灰度級z具有規定 PDF 的影像，

2.2.4 、連續灰度影像的直方圖規定化小結

通過前面的原理及案例介紹，我們看到，連續灰度值的非數字影像的直方圖規定化是很容易理解的，并且均衡輸人影像的中間一步可以跳過；我們需要的是得到將r映射為s的變換函式T（r），然后，這兩步可以合并為從r到z的一步變換，即由：

就可以求出T（r），

但實際中，比較困難的是尋找T（r）和G的反函式形式，老猿認為由于輸出影像的PDF是指定的，要找T（r）和G的反函式其實最困難的是一副輸入影像的pr（r）怎么求得，因為輸入影像只有各個影像元素的灰度值，

2.3、數字影像的直方圖匹配或規定化

2.3.1、相關原理

連續灰度非數字影像求其映射到結果影像的映射函式比較困難（老猿認為主要是輸入影像的PDF函式不好求），但幸運的是，計算機處理的是數字影像，因此相關解決方案可以大大簡化，

在《數字影像處理：直方圖均衡（Histogram Equalization）的原理及處理介紹》中介紹了數字影像直方圖均衡變換的公式：

它是式（3.3-10）的離散形式，其中MN是影像像素的行數M與列數N的乘積，L為該表示影像灰度所在數字區間的上限，即影像中可能的灰度級數（intensity levels），如8bit影像為256，nj表示灰度值=j的像素個數，sk為輸入影像灰度值為rk經過變換后映射到目標影像的灰度值，T（rk）表示將輸入灰度值rk映射到sk的變換函式，

假設有一幅與輸入影像內容相同灰度值不同的數字影像，假設其所有灰度值為Z0–Zq，每個灰度值zi出現的概率為pz（zi），那么對其進行直方圖均衡對應的公式為式（3.3-11）對應的離散形式：

由于G(Zq)=Sk，因此有：

也就是說，對于灰度值為Sk的輸入影像，對任意給定了輸出影像的灰度值Zq，能找到一個從輸入影像的灰度值映射到輸出影像灰度值的變換函式，從而實作將輸入影像映射到指定灰度值的輸出影像，注意這里的Sk對應的灰度值為輸入影像灰度為rk的灰度經直方圖均衡變換后的灰度值，

實踐中使用該方法處理數字影像時，由于數字影像灰度值的有限性，我們可以通過式（3.3-14）計算出指定灰度值的輸出影像的直方圖均衡后的近似灰度值（因為數字影像灰度值為整數，因此對式（3.3-14）的計算結果要進行四舍五入），從而我們就構建了一個Zq到Sk的映射表，這樣無需求出G函式的反函式，直接通過映射表的Sk的取值查到Zq的取值，例如表中第64行記錄的Sk值最接近給定的Sk值，則Zq的q值=63（注意是從0開始計數），故Z63是G(Zq)=Sk最近似解，

注意：使用式（3.3-14）計算出來的映射表，其資料源值為Z0-Zn，即對應灰度值范圍為[0,L-1]的所有升序排列的整數，即Z0=0、Z1=1、…、n=l-1，其結果值為Sk對應的值，由于直方圖均衡變換四舍五入取整后有可能Sk值組成的集合不包含[0,L-1]內所有整數，

2.3.2、數字影像直方圖規定化處理步驟

假設Sk是直方圖均衡后影像的灰度值，可以總結直方圖規定化程序如下：

1. 計算給定影像的直方圖pr（r），并用它尋找式(3.3-13)的直方圖均衡變換，把變換后影像的灰度sk四舍五入為范圍[0，L-1]內的整數；
2. 用式(3.3-14)對q=0、1、2、…、L-1計算變換函式G的所有值，其中pz(zi)是規定的直方圖的值，把G的值四舍五入為范圍[0，L-1]內的整數，并將G的輸入和結果值存盤在一個表中，注意式(3.3-14)是對規定的直方圖進行直方圖均衡變換，其結果為均衡化后直方圖的灰度值，其值與步驟1中的sk一致，而輸入是zi，zi的灰度值為i；
3. 對每一個值sk，k=0、1、2、…、L-1，使用步驟2存盤的G值尋找相應的zq值，以使G(zq)最接近sk，并存盤這些從s到z的映射，當滿足給定sk的zq值多于一個時(即映射不唯一時)，按慣例選擇最小的值；
4. 首先對輸入影像進行均衡，然后使用步驟 3 找到的映射把該影像中的每個均衡后的像素值sk映射為直方圖規定化后的影像中的相應zq值，形成直方圖規定化后的影像，正如連續情況那樣，均衡輸入影像的中間步驟是概念上的，它可以用合并兩個變換函式T和G的反函式跳過這一步，

如早些時候提到的那樣，對于滿足條件(a’)嚴格單調和(b)函式值屬于[0,L-1]范圍的G的反函式，G必須是嚴格單調的，根據式(3.3-14)，它意味著規定直方圖的任何pz(zi)值都不能為零，當作業在離散數值的情況時，該條件可能不滿足的事實并不是一個嚴重的實作問題，采用步驟3中取最小值的方式就可以解決，

注意：步驟4前面幾句的翻譯有問題，看起來莫名其妙，其對應英文原文如下：
Form the histogram-specified image by first histogram-equalizing the input image and then mapping every equalized pixel value, sk,of this image to the corresponding value zq in the histogram-specified image using the mappings found in step 3.

老猿將其翻譯在付費文章《數字影像處理：直方圖匹配或規定化Histogram Matching (Specification)原理及感悟》中進行了修改，上述處理步驟實際上還是不怎么好理解，老猿也在該文中按自己的理解進行了介紹，

2.3.3、數字影像的直方圖規定化案例

假設一幅大小為64x64像素(MN=4096)的3位元影像(L=8)的灰度分布如表3.1所示，其中灰度級是范圍[0，L-1]=[0，7]中的整數，其灰度分布和直方圖值如下：

現要求變換該直方圖，以便使其具有表3.2第2列規定的值：

第一步，利用式(3.3-13)計算該輸入影像均衡后的sk值，四舍五入后可以得到輸入影像和均衡化影像的映射關系，sk中的k即為輸入影像的對應灰度值：

第二步，利用式(3.3-14)計算表3.2中第2列的值計算規定直方圖與均衡化后影像灰度值sk的映射關系，四舍五入后可以得到：

即得到sk（G(zq)）和zq的映射關系：

使用步驟3中的方式查找每個sk對應的最小的zq值：

這兒有個問題，為什么sk為2的值被替換為了3？老猿才開始以為是印刷問題，看了英文原版也是3，仔細思考了下，終于理解了，相關介紹老猿在《數字影像處理：直方圖匹配或規定化Histogram Matching (Specification)原理及感悟》中進行了介紹，

得到表3.4以后，我們使用表中的映射把直方圖均衡后的影像中的每個像素映射為新創建的直方圖規定化影像中的相應像素，這樣得到的結果直方圖的值列在表3.2的第三列，由于sk值為2的值被替換為了3，導致最終結果并不完全與規定直方圖匹配，但達到了將輸入影像灰度朝指定灰度高度近似的效果，

上述處理步驟中，得到直方圖均衡后影像的中間步驟有助于解釋該匹配程序，但這個步驟不是必須的，只需要得到r到s的映射和s到z的映射，然后用這些映射把原始像素直接映射到直方圖規定化后的影像像素即可，

2.3.4、直方圖均衡與直方圖匹配案例對比

上圖3.23(a)顯示了火星衛星(Phobos)影像，該影像由 NASA 的火星表面探測器拍攝，圖3.23(b)顯示了圖3.23(a)的直方圖，影像的大部分是暗色區域，產生了像素集中于灰度級暗端的直方圖的特點，乍看之下，會得到直方圖均衡化是增強該影像的一個好辦法的結論，以便使暗區域的細節更清楚，下面的討論將證明事實并非如此，

上圖3.24(a)顯示了由圖3.23(b)所示直方圖得到的直方圖均衡變換[式(3.3-8)或式(3.3-13)]，這個變換函式最相關的特征是它從灰度0快速上升到接近190，這是因為輸人直方圖中有大量的像素灰度接近0，當這種變換應用于輸入影像的灰度，以獲得直方圖的均衡效果時，最終結果是把非常窄的暗像素區間映射到輸出影像灰度級的高端，因為輸入影像的大量像素在此暗像素區間有精確的灰度值，所以我們希望結果是具有明亮、“沖淡”（washed-out）外觀的影像，如圖3.24(b)所示，情況確實如此，該影像的直方圖示于圖3.24?，注意，所有灰度級基本上都偏向了灰度范圍的上半部，

由于圖3.24(a)中變換函式的問題是由原影像中大量灰度接近于0的像素引起的，一種合理的方法是修改影像的直方圖，以便使其不再有這一屬性，

圖3.25(a)顯示了一個人為規定的函式，該函式保持了原有直方圖（圖3.23（b））的基本形狀，但在影像灰度暗區的灰度值卻有更平滑的過渡，將這個函式取樣為256個等間隔
的離散值，就生成了所需的規定的直方圖，

使用式(3.3-14)從這一直方圖得到的變換函式G(z)，該變換在圖3.25(b)中標注為變換(1)，類似地，由式(3.3-16)(使用前面討論的逐步程序得到)匯出的反變換函式G^-1(s)在圖3.25(b)中標注為變換(2)，

圖3.25?中增強后的影像是把變換(2)用于圖3.24(b)中直方圖均衡化后的影像的像素得到的，比較這兩幅影像可以看出，在使用直方圖均衡得到的結果的基礎上，直方圖規定化影像的改進是很明顯的，

注意，原始直方圖的細小改變會使影像的外觀得到顯著改進，圖3.25(d)顯示了圖3.25?的直方圖，該直方圖最容易辨別的特征是其低端已向右移向灰度級的較亮區域(但并不過分)，這正是我們所希望的，

三、小結

本文介紹了連續影像和數字影像直方圖匹配（直方圖規定化）的原理、處理程序，并提供了案例進行了講解，可以看到，直方圖均衡處理是直方圖匹配的一個重要橋梁，

到目前為止，雖然我們對直方圖規定化已經很清楚了，但在結束這一節之前，我們要強調直方圖規定化在大多數時候都是試錯（trial-and-error）程序，這正像前邊的例子中我們所做的那樣，一個可用的指導路線是用手邊的問題學習，但是，有時可能會遇到這樣的情況，即定義“平均”直方圖的外觀，并把它作為規定的直方圖，在這種情況下，直方圖規定化就變成了一種直接的處理，然而，一般來說，并沒有規定直方圖的規則，對于任何一個給定的增強任務都必須借助于實際分析，

最后，對于直方圖規定化這個翻譯個人覺得很low，個人感覺叫指定直方圖更好理解，但Histogram Specification怎么翻譯確實不太好把握，因此還是叫Histogram Matching直方圖匹配最好，

更多直方圖處理相關知識請參考《《數字影像處理》第三章學習總結感悟2：直方圖處理》，

更多影像處理請參考專欄《OpenCV-Python圖形影像處理》及《影像處理基礎知識》的介紹，

對于缺乏Python基礎的同仁，可以通過老猿的免費專欄《 專欄：Python基礎教程目錄》從零開始學習Python，

如對文章內容存在疑問，可在博客評論區留言，或關注：老猿Python 微信公號發訊息咨詢，可通過掃描博客左邊的二維碼加微信公眾號，

寫博不易，敬請支持：

如果閱讀本文于您有所獲，敬請點贊、評論、收藏，謝謝大家的支持！

關于老猿的付費專欄

1. 付費專欄《https://blog.csdn.net/laoyuanpython/category_9607725.html 使用PyQt開發圖形界面Python應用》專門介紹基于Python的PyQt圖形界面開發基礎教程，對應文章目錄為《 https://blog.csdn.net/LaoYuanPython/article/details/107580932 使用PyQt開發圖形界面Python應用專欄目錄》；
2. 付費專欄《https://blog.csdn.net/laoyuanpython/category_10232926.html moviepy音視頻開發專欄 )詳細介紹moviepy音視頻剪輯合成處理的類相關方法及使用相關方法進行相關剪輯合成場景的處理，對應文章目錄為《https://blog.csdn.net/LaoYuanPython/article/details/107574583 moviepy音視頻開發專欄文章目錄》；
3. 付費專欄《https://blog.csdn.net/laoyuanpython/category_10581071.html OpenCV-Python初學者疑難問題集》為《https://blog.csdn.net/laoyuanpython/category_9979286.html OpenCV-Python圖形影像處理 》的伴生專欄，是筆者對OpenCV-Python圖形影像處理學習中遇到的一些問題個人感悟的整合，相關資料基本上都是老猿反復研究的成果，有助于OpenCV-Python初學者比較深入地理解OpenCV，對應文章目錄為《https://blog.csdn.net/LaoYuanPython/article/details/109713407 OpenCV-Python初學者疑難問題集專欄目錄 》
4. 付費專欄《https://blog.csdn.net/laoyuanpython/category_10762553.html Python爬蟲入門 》站在一個互聯網前端開發小白的角度介紹爬蟲開發應知應會內容，包括爬蟲入門的基礎知識，以及爬取CSDN文章資訊、博主資訊、給文章點贊、評論等實戰內容，

前兩個專欄都適合有一定Python基礎但無相關知識的小白讀者學習，第三個專欄請大家結合《https://blog.csdn.net/laoyuanpython/category_9979286.html OpenCV-Python圖形影像處理 》的學習使用，

對于缺乏Python基礎的同仁，可以通過老猿的免費專欄《https://blog.csdn.net/laoyuanpython/category_9831699.html 專欄：Python基礎教程目錄）從零開始學習Python，

如果有興趣也愿意支持老猿的讀者，歡迎購買付費專欄，

老猿Python，跟老猿學Python！

? ? 前往老猿Python博文目錄 https://blog.csdn.net/LaoYuanPython ?