基于直方圖的圖像增強(qiáng)算法（HE、CLAHE、Retinex）之（二）

　　作為圖像增強(qiáng)算法系列的第二篇文章，下面我們將要介紹功能強(qiáng)大、用途廣泛、影響深遠(yuǎn)的對比度有限的自適應(yīng)直方圖均衡(CLAHE，Contrast Limited Adaptive Histogram Equalization)算法。盡管最初它僅僅是被當(dāng)作一種圖像增強(qiáng)算法被提出，但是現(xiàn)今在圖像去霧、低照度圖像增強(qiáng)，水下圖像效果調(diào)節(jié)、以及數(shù)碼照片改善等方面都有應(yīng)用。這個(gè)算法的算法原理看似簡單，但是實(shí)現(xiàn)起來卻并不那么容易。我們將結(jié)合相應(yīng)的Matlab代碼來對其進(jìn)行解釋。希望你在閱讀本文之前對樸素的直方圖均衡算法有所了解，相關(guān)內(nèi)容可以參見本系列的前一篇文章：基于直方圖的圖像增強(qiáng)算法(HE、CLAHE、Retinex)之(一)

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

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

　　效果圖A 效果圖B

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

[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)

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

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

　　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);

　　上述程序之結(jié)果效果圖A所示。

　　下面程序?qū)⒃瓐D像的色彩空間轉(zhuǎn)換到LAB空間之后再對L通道進(jìn)行處理。

　　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);

　　上述程序所得之結(jié)果如圖B所示。

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

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

　　clc;

　　clear all;

　　Img = rgb2gray(imread('space.jpg'));

　　[h,w] = size(Img);

　　minV = double(min(min(Img)));

　　maxV = double(max(max(Img)));

　　imshow(Img);

　　圖像的初始狀態(tài)顯示如下。此外該圖的 Height = 395，Width = 590，灰度最大值為255，最小值為8。

　　我們希望把原圖像水平方向分成8份，把垂直方向分成4份，即原圖將被劃分成4 × 8 = 32個(gè)SubImage。然后可以算得每個(gè)塊(tile)的height = 99，width = 74。注意，由于原圖的長、寬不太可能剛好可被整除，所以我在這里的處理方式是建立一個(gè)稍微大一點(diǎn)的圖像，它的寬和長都被補(bǔ)上了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;

　　對長和寬進(jìn)行填補(bǔ)之后，對新圖像的一些必要信息進(jìn)行更新。

　　new_w = w + deltax;

　　new_h = h + deltay;

　　NrPixels = WSize * WSize;

　　然后指定圖像中直方圖橫坐標(biāo)上取值的計(jì)數(shù)(也就指定了統(tǒng)計(jì)直方圖上橫軸數(shù)值的間隔或計(jì)數(shù)的精度)，對于色彩比較豐富的圖像，我們一般都要求這個(gè)值應(yīng)該大于128。

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

　　NrBins = 256;

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

　　LUT = zeros(maxV+1,1);

　　for i=minV:maxV

　　LUT(i+1) = fix(i - minV);%i+1

　　end

　　Bin = zeros(new_h, new_w);

　　for m = 1 : new_h

　　for n = 1 : new_w

　　Bin(m,n) = 1 + LUT(tmpImg(m,n) + 1);

　　end

　　end

　　Hist = zeros(NrY, NrX, 256);

　　for i=1:NrY

　　for j=1:NrX

　　tmp = 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

　　end

　　Hist = circshift(Hist,[0, 0, -1]);

　　注意：按通常的理解，上面這一步我們應(yīng)該建立的直方圖(集合)應(yīng)該是一個(gè)4×8=32個(gè)長度為256的向量(你當(dāng)然也可以這么做)。但由于涉及到后續(xù)的一些處理方式，我這里是生成了一個(gè)長度為256的4×8矩陣。Index = 1的矩陣其實(shí)相當(dāng)于是整張圖像各個(gè)tile上灰度值=0的像素個(gè)數(shù)計(jì)數(shù)。例如，我們所得的Hist(:, :, 18)如下。這就表明圖像中最左上角的那個(gè)tile里面灰度值=17的像素有零個(gè)。同理，它右邊的一個(gè)tile則有46個(gè)灰度值=17的像素。

　　Hist(:,:,18) =

　　0 46 218 50 14 55 15 7

　　0 0 21 18 114 15 74 73

　　0 1 0 0 2 67 124 82

　　0 0 0 0 0 1 9 2

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

　　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);

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

　　[plain] view plain copytmp_hist = reshape(Hist(1,1,:), 1, 256);

　　plot(tmp_hist)

　　輸出結(jié)果下圖中的左圖所示。

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

　　繼續(xù)回到程序?qū)崿F(xiàn)上的討論。最后，也是最關(guān)鍵的步驟，我們需要對結(jié)果進(jìn)程插值處理。這也是Zuiderveld設(shè)計(jì)的算法中最復(fù)雜的部分。

　　yI = 1;

　　for i = 1:NrY+1

　　if i == 1

　　subY = floor(HSize/2);

　　yU = 1;

　　yB = 1;

　　elseif i == NrY+1

　　subY = floor(HSize/2);

　　yU = NrY;

　　yB = NrY;

　　else

　　subY = HSize;

　　yU = i - 1;

　　yB = i;

　　end

　　xI = 1;

　　for j = 1:NrX+1

　　if j == 1

　　subX = floor(WSize/2);

　　xL = 1;

　　xR = 1;

　　elseif j == NrX+1

　　subX = floor(WSize/2);

　　xL = NrX;

　　xR = NrX;

　　else

　　subX = WSize;

　　xL = j - 1;

　　xR = j;

　　end

　　UL = 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 - 1

　　inverseI = subY - i;

　　for j = 0:subX - 1

　　inverseJ = 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;

　　end

　　end

　　%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

　　output(yI:yI+subY-1, xI:xI+subX-1) = sImage;

　　xI = xI + subX;

　　end

　　yI = yI + subY;

　　end

　　這個(gè)地方，作者原文中已經(jīng)講得比較清楚了，我感覺我也沒有必要狗尾續(xù)貂，班門弄斧了。下面截作者原文中的一段描述，足以說明問題。

　　最后來看看我們處理的效果如何(當(dāng)然，這里還需要把之前我們填補(bǔ)的部分裁掉)。

　　output = output(1:h, 1:w);

　　figure, imshow(output, []);

　　來看看結(jié)果吧~可以對比一下之前的灰度圖，不難發(fā)現(xiàn)，圖像質(zhì)量已有大幅改善。