基於直方圖的圖像增強算法(HE、CLAHE、Retinex)之(二)

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

本文引用地址：http://www.eepw.com.cn/article/201702/344364.htm

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

　　下面是利用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: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]);

　　注意：按通常的理解，上面這一步我們應該建立的直方圖(集合)應該是一個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 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

　　然後來對直方圖進行裁剪。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為例，它的原直方圖分布可以用下面代碼來繪出：

　　[plain] view plain copytmp_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+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

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

　　最後來看看我們處理的效果如何(當然，這裡還需要把之前我們填補的部分裁掉)。

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

　　figure, imshow(output, []);

　　來看看結果吧~可以對比一下之前的灰度圖，不難發現，圖像質量已有大幅改善。