直方圖是圖像色彩統計特征的抽象表述。基于直方圖可以實現很多有趣的算法。例如，圖像增強中利用直方圖來調整圖像的對比度、有人利用直方圖來進行大規模無損數據隱藏、還有人利用梯度直方圖HOG來構建圖像特征進而實現目標檢測。本節我們就來討論重要的直方圖均衡化算法，說它重要是因為以此為基礎后續又衍生出了許多實用而有趣的算法。

Histogram equalization

如果一幅圖像的像素灰度值在一個過于有限的范圍內聚集，那么圖像的程序效果即會很糟糕，直接觀感就是對比度很弱。下圖來自維基百科，第一幅圖的直方圖分布非常不均衡。如果把直方圖均勻地延展到整個分布域內，則圖像的效果顯得好了很多。

Matlab中提供了現成的函數“histeq()”來實現圖像的直方圖均衡。但為了演示說明算法的原理，下面我將在Matlab中自行編碼實現圖像的直方圖均衡。通過代碼來演示這個算法顯然更加直觀，更加易懂。當然，其實我還不得不感嘆，如果僅僅是作為圖像算法研究之用，Matlab確實非常好用。

首先讀入圖像，并將其轉化為灰度圖。然后提取圖像的長和寬。

image = imread('Unequalized_Hawkes_Bay_NZ.jpg'); Img = rgb2gray(image); [height,width]=size(image);
然后繪制一下原始圖像的直方圖。

[counts1, x] = imhist(Img,256); counts2 = counts1/height/width; stem(x, counts2);

統計每個灰度的像素值累計數目。

NumPixel = zeros(1,256);%統計各灰度數目，共256個灰度級 for i = 1:heightfor j = 1: width%對應灰度值像素點數量增加一%因為NumPixel的下標是從1開始，但是圖像像素的取值范圍是0~255，所以用NumPixel(Img(i,j) + 1)NumPixel(Img(i,j) + 1) = NumPixel(Img(i,j) + 1) + 1;end end
然后將頻數值算為頻率

ProbPixel = zeros(1,256); for i = 1:256ProbPixel(i) = NumPixel(i) / (height * width * 1.0); end
再用函數cumsum來計算cdf，并將頻率（取值范圍是0.0~1.0）映射到0~255的無符號整數。

CumuPixel = cumsum(ProbPixel); CumuPixel = uint8(255 .* CumuPixel + 0.5);

直方圖均衡。賦值語句右端，Img(i,j)被用來作為CumuPixel的索引。比如Img(i,j)?= 120，則從CumuPixel中取出第120個值作為Img(i,j)?的新像素值。

for i = 1:heightfor j = 1: widthImg(i,j) = CumuPixel(Img(i,j));end end

最后顯示新圖像的直方圖。

imshow(Img); [counts1, x] = imhist(Img,256); counts2 = counts1/height/width; stem(x, counts2);

當然，上述討論的是灰度圖像的直方圖均衡。對于彩色圖像而言，你可能會想到分別對R、G、B三個分量來做處理，這也確實是一種方法。但有些時候，這樣做很有可能導致結果圖像色彩失真。因此有人建議將RGB空間轉換為HSV之后，對V分量進行直方圖均衡處理以保存圖像色彩不失真。下面我們來做一些對比實驗。待處理圖像是標準的圖像處理測試用圖couple圖，如下所示。

首先，我們分別處理R、G、B三個分量，為了簡便我們直接使用matlab中的函數histeq()。

a = imread('couple.tiff'); R = a(:,:,1); G = a(:,:,2); B = a(:,:,3);R = histeq(R, 256); G = histeq(G, 256); B = histeq(B, 256);a(:,:,1) = R; a(:,:,2) = G; a(:,:,3) = B; imshow(a)
下面的代碼使用了另外一種方式，即將色彩空間轉換到HSV后，對V通道進行處理。由于代碼基本與前面介紹的一致，這里我們不再做過多解釋了。

Img = imread('couple.tiff'); hsvImg = rgb2hsv(Img); V=hsvImg(:,:,3); [height,width]=size(V);V = uint8(V*255); NumPixel = zeros(1,256); for i = 1:heightfor j = 1: widthNumPixel(V(i,j) + 1) = NumPixel(V(i,j) + 1) + 1;end endProbPixel = zeros(1,256); for i = 1:256ProbPixel(i) = NumPixel(i) / (height * width * 1.0); endCumuPixel = cumsum(ProbPixel); CumuPixel = uint8(255 .* CumuPixel + 0.5);for i = 1:heightfor j = 1: widthV(i,j) = CumuPixel(V(i,j));end endV = im2double(V); hsvImg(:,:,3) = V; outputImg = hsv2rgb(hsvImg); imshow(outputImg);

最后，來對比一下不同方法對彩色圖像的處理效果。下面的左圖是采用R、G、B三分量分別處理得到的結果。右圖是對HSV空間下V通道處理之結果。顯然，右圖的效果更理想，而左圖則出現了一定的色彩失真。事實上，對彩色圖像進行直方圖均衡是圖像處理研究領域一個看似簡單，但是一直有人在研究的話題。我們所說的對HSV空間中V分量進行處理的方法也是比較基本的策略。很多相關的研究文章都提出了更進一步的、適應性更強的彩色圖像直方圖均衡化算法。有興趣的讀者可以參閱相關文獻以了解更多。

分別處理R、G、B三個分量之結果??? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ?? 轉換到HSV空間后處理V分量

自適應直方圖均衡（CLAHE，Contrast Limited Adaptive Histogram Equalization）算法

先來看一下待處理的圖像效果：

下面是利用CLAHE算法處理之后得到的兩個效果（后面我們還會具體介紹我們所使用的策略）：

?? ?

???????????????????????????? 效果圖A???????????????????????????????????????????????????????????????????????? 效果圖B

對于一幅圖像而言，它不同區域的對比度可能差別很大。可能有些地方很明亮，而有些地方又很暗淡。如果采用單一的直方圖來對其進行調整顯然并不是最好的選擇。于是人們基于分塊處理的思想提出了自適應的直方圖均衡算法AHE。維基百科上說的也比較明白：AHE improves on this by transforming each pixel with a transformation function derived from a neighbourhood region. 但是這種方法有時候又會將一些噪聲放大，這是我們所不希望看到的。于是荷蘭烏得勒支大學的Zuiderveld教授又引入了CLAHE，利用一個對比度閾值來去除噪聲的影響。特別地，為了提升計算速度以及去除分塊處理所導致的塊邊緣過渡不平衡效應，他又建議采用雙線性插值的方法。關于算法的介紹和描述，下面這兩個資源已經講得比較清楚。

[1] https://en.wikipedia.org/wiki/Adaptive_histogram_equalization#Contrast_Limited_AHE

[2] K. Zuiderveld:?Contrast Limited Adaptive Histogram Equalization. In: P. Heckbert:?Graphics Gems IV, Academic Press 1994 （http://www.docin.com/p-119164091.html）

事實上，盡管這個算法原理，然而它實現起來卻仍然有很多障礙。但在此之前，筆者還需說明的是，Matlab中已經集成了實現CLAHE的函數adapthisteq()，如果你僅僅需要一個結果，其實直接使用這個函數就是最好的選擇。我給出一些示例代碼用以生成前面給出之效果。函數adapthisteq()只能用來處理灰度圖，如果要處理彩色圖像，則需要結合自己編寫的代碼來完成。上一篇文章介紹了對彩色圖像進行直方圖均衡的兩種主要策略：一種是對R、G、B三個通道分別進行處理；另外一種是轉換到另外一個色彩空間中再進行處理，例如HSV（轉換后只需對V通道進行處理即可）。

首先，我們給出對R、G、B三個通道分別使用adapthisteq()函數進行處理的示例代碼：

img = imread('space.jpg'); rimg = img(:,:,1); gimg = img(:,:,2); bimg = img(:,:,3); resultr = adapthisteq(rimg); resultg = adapthisteq(gimg); resultb = adapthisteq(bimg); result = cat(3, resultr, resultg, resultb); imshow(result); 上述程序之結果效果圖A所示。

下面程序將原圖像的色彩空間轉換到LAB空間之后再對L通道進行處理。

clear; img = imread('space.jpg'); cform2lab = makecform('srgb2lab'); LAB = applycform(img, cform2lab); L = LAB(:,:,1); LAB(:,:,1) = adapthisteq(L); cform2srgb = makecform('lab2srgb'); J = applycform(LAB, cform2srgb); imshow(J);
上述程序所得之結果如圖B所示。

如果你希望把這個算法進一步提升和推廣，利用用于圖像去霧、低照度圖像改善和水下圖像處理，那么僅僅知其然是顯然不夠的，你還必須知其所以然。希望我下面一步一步實現的代碼能夠幫你解開這方面的困惑。鑒于前面所列之文獻已經給出了比較詳細的算法描述，下面將不再重復這部分內容，轉而采用Matlab代碼來對其中的一些細節進行演示。

首先來從灰度圖的CLAHE處理開始我們的討論。為此清理一下Matlab的環境。然后，讀入一張圖片（并將其轉化灰度圖），獲取圖片的長、寬、像素灰度的最大值、最小值等信息。

clc; clear all; Img = rgb2gray(imread('space.jpg')); [h,w] = size(Img); minV = double(min(min(Img))); maxV = double(max(max(Img))); imshow(Img);
圖像的初始狀態顯示如下。此外該圖的 Height = 395，Width = 590，灰度最大值為255，最小值為8。

我們希望把原圖像水平方向分成8份，把垂直方向分成4份，即原圖將被劃分成4 × 8 = 32個SubImage。然后可以算得每個塊（tile）的height = 99，width = 74。注意，由于原圖的長、寬不太可能剛好可被整除，所以我在這里的處理方式是建立一個稍微大一點的圖像，它的寬和長都被補上了deltax和deltay，以保證長、寬都能被整除。

NrX = 8; NrY = 4; HSize = ceil(h/NrY); WSize = ceil(w/NrX);deltay = NrY*HSize - h; deltax = NrX*WSize - w;tmpImg = zeros(h+deltay,w+deltax); tmpImg(1:h,1:w) = Img;

對長和寬進行填補之后，對新圖像的一些必要信息進行更新。

new_w = w + deltax; new_h = h + deltay;NrPixels = WSize * WSize;
然后指定圖像中直方圖橫坐標上取值的計數（也就指定了統計直方圖上橫軸數值的間隔或計數的精度），對于色彩比較豐富的圖像，我們一般都要求這個值應該大于128。

% NrBins - Number of greybins for histogram ("dynamic range") NrBins = 256;

然后用原圖的灰度取值范圍重新映射了一張Look-Up Table（當然你也可以直接使用0~255這個范圍，這取決你后續建立直方圖的具體方法），并以此為基礎為每個圖像塊（tile）建立直方圖。

LUT = zeros(maxV+1,1);for i=minV:maxVLUT(i+1) = fix(i - minV);%i+1 endBin = zeros(new_h, new_w); for m = 1 : new_hfor n = 1 : new_wBin(m,n) = 1 + LUT(tmpImg(m,n) + 1);end endHist = zeros(NrY, NrX, 256); for i=1:NrYfor j=1:NrXtmp = uint8(Bin(1+(i-1)*HSize:i*HSize, 1+(j-1)*WSize:j*WSize));%tmp = tmpImg(1+(i-1)*HSize:i*HSize,1+(j-1)*WSize:j*WSize);[Hist(i, j, :), x] = imhist(tmp, 256);end endHist = circshift(Hist,[0, 0, -1]);
注意：按通常的理解，上面這一步我們應該建立的直方圖（集合）應該是一個4×8=32個長度為256的向量（你當然也可以這么做）。但由于涉及到后續的一些處理方式，我這里是生成了一個長度為256的4×8矩陣。Index = 1的矩陣其實相當于是整張圖像各個tile上灰度值=0的像素個數計數。例如，我們所得的Hist(:, :, 18)如下。這就表明圖像中最左上角的那個tile里面灰度值=17的像素有零個。同理，它右邊的一個tile則有46個灰度值=17的像素。

Hist(:,:,18) =0 46 218 50 14 55 15 70 0 21 18 114 15 74 730 1 0 0 2 67 124 820 0 0 0 0 1 9 2

然后來對直方圖進行裁剪。Matlab中內置的函數adapthisteq()中cliplimit參數的取值范圍是0~1。這里我們所寫的方法則要求該值>1。當然這完全取決你算法實現的策略，它們本質上并沒有差異。然后我們將得到新的（裁剪后的）映射直方圖。

ClipLimit = 2.5; ClipLimit = max(1,ClipLimit * HSize * WSize/NrBins);Hist = clipHistogram(Hist,NrBins,ClipLimit,NrY,NrX); Map=mapHistogram(Hist, minV, maxV, NrBins, NrPixels, NrY, NrX);

因為這里沒有具體給出clipHistogram函數的實現，所以此處我希望插入一部分內容來解釋一下我的實現策略（也就是說，在實際程序中并不需要包含這部分）。我們以圖像最左上角的一個tile為例，它的原直方圖分布可以用下面代碼來繪出：

tmp_hist = reshape(Hist(1,1,:), 1, 256); plot(tmp_hist)

輸出結果下圖中的左圖所示。

如果我們給ClipLimit賦初值為2.5，則經過語句ClipLimit = max(1,ClipLimit * HSize * WSize/NrBins);計算之后，ClipLimit將變成71.54。然后我們再用上述代碼繪制新的直方圖，其結果將如上圖中的右圖所示。顯然，圖中大于71.54的部分被裁剪掉了，然后又平均分配給整張直方圖，所以你會發現整張圖都被提升了。這就是我們這里進行直方圖裁剪所使用的策略。但是再次強調，matlab中的內置函數adapthisteq()僅僅是將這個參數進行了歸一化，這與我們所使用的方法并沒有本質上的區別。

繼續回到程序實現上的討論。最后，也是最關鍵的步驟，我們需要對結果進程插值處理。這也是Zuiderveld設計的算法中最復雜的部分。

yI = 1; for i = 1:NrY+1if i == 1subY = floor(HSize/2);yU = 1;yB = 1;elseif i == NrY+1subY = floor(HSize/2);yU = NrY;yB = NrY;elsesubY = HSize;yU = i - 1;yB = i;endxI = 1;for j = 1:NrX+1if j == 1subX = floor(WSize/2);xL = 1;xR = 1;elseif j == NrX+1subX = floor(WSize/2);xL = NrX;xR = NrX;elsesubX = WSize;xL = j - 1;xR = j;endUL = Map(yU,xL,:);UR = Map(yU,xR,:);BL = Map(yB,xL,:);BR = Map(yB,xR,:);subImage = Bin(yI:yI+subY-1,xI:xI+subX-1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sImage = zeros(size(subImage));num = subY * subX;for i = 0:subY - 1inverseI = subY - i;for j = 0:subX - 1inverseJ = subX - j;val = subImage(i+1,j+1);sImage(i+1, j+1) = (inverseI*(inverseJ*UL(val) + j*UR(val)) ...+ i*(inverseJ*BL(val) + j*BR(val)))/num;endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%output(yI:yI+subY-1, xI:xI+subX-1) = sImage;xI = xI + subX;endyI = yI + subY; end
這個地方，作者原文中已經講得比較清楚了，我感覺我也沒有必要狗尾續貂，班門弄斧了。下面截作者原文中的一段描述，足以說明問題。

最后來看看我們處理的效果如何（當然，這里還需要把之前我們填補的部分裁掉）。

output = output(1:h, 1:w); figure, imshow(output, []);
來看看結果吧~可以對比一下之前的灰度圖，不難發現，圖像質量已有大幅改善。

總結

以上是生活随笔為你收集整理的基于直方图的图像增强算法（HE、CLAHE、Retinex）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。