對比度受限的自適應直方圖均衡化(CLAHE)

直方圖均衡化（HE）是一種很常用的直方圖類方法，基本思想是通過圖像的灰度分布直方圖確定一條映射曲線，用來對圖像進行灰度變換，以達到提高圖像 對比度的目的。該映射曲線其實就是圖像的累計分布直方圖（CDF）（嚴格來說是呈正比例關系）。然而HE是對圖像全局進行調整的方法，不能有效地提高局部 對比度，而且某些場合效果會非常差。如：

上述原圖和HE結果圖的直方圖分別為：

因為從原圖的直方圖中求取的映射函數（CDF）形狀為：

將它作用於原圖像會導致直方圖被整體右移，沒有充分利用整個灰度動態范圍。

為了提高圖像的局部對比度，有人提出將圖像分成若干子塊，對子塊進行HE處理，這便是AHE（自適應直方圖均衡化），使用AHE處理上圖得到：

 

結果直方圖：

可 以看出結果圖像的灰度較好地分布在了全部動態范圍上。從結果圖像上也可以看出，局部對比度的確得到了提高，視覺效果要優於HE。但是仍然有個問題：AHE 對局部對比度提高過大，導致圖像失真。看看背景區，本來的黑色背景現在已經變成白色了，原因是因為背景區中的局部子塊統計得到的直方圖在0灰度處幅值太高 （實際上全黑子圖基本上就集中在0灰度處），這樣導致映射曲線斜率過高，將所有灰度值都映射到整個灰度軸的右側，所以結果圖中背景偏白（另外局部對比度過 高還會放大圖像中的噪聲，不過上圖並沒有體現這一點）。

為了解決這個問題，我們必須對局部對比度進行限制，這就是我們今天的主題：CLAHE。

從HE中我們知道，映射曲線T與CDF關系為：（M為最高灰度值，N為像素個數）

限制對比度，其實就是限制CDF的斜率，又因累計分布直方圖CDF是灰度直方圖Hist的積分：

反過來：

也就是說限制CDF的斜率就相當於限制Hist的幅度。

因此我們需要對子塊中統計得到的直方圖進行裁剪，使其幅值低於某個上限，當然裁剪掉的部分又不能扔掉，我們還需要將這部分裁剪值均勻地分布在整個灰度區間上，以保證直方圖總面積不變，如下圖：

可以看到，這時直方圖又會整體上升了一個高度，貌似會超過我們設置的上限。其實在具體實現的時候有很多解決方法，你可以多重復幾次裁剪過程，使得上升的部分變得微不足道，或是用另一種常用的方法：

設 裁剪值為ClipLimit，求直方圖中高於該值的部分的和totalExcess，此時假設將totalExcess均分給所有灰度級，  求出這樣導致的直方圖整體上升的高度L=totalExcess/N，以upper= ClipLimit-L為界限對直方圖進行如下處理：

(1)若幅值高於ClipLimit，直接置為ClipLimit；

(2)若幅值處於Upper和ClipLimit之間，將其填補至ClipLimit；

(3)若幅值低於Upper，直接填補L個像素點；

經過上述操作，用來填補的像素點個數通常會略小於totalExcess，也就是還有一些剩余的像素點沒分出去，這個剩余來自於(1)(2)兩處。這時我們可以再把這些點均勻地分給那些目前幅值仍然小於ClipLimit的灰度值。

這里給出一段代碼：（摘自Matlab的adapthisteq.m），描述的就是上述過程：

% total number of pixels overflowing clip limit in each bin
totalExcess = sum(max(imgHist - clipLimit,0));  
 
% clip the histogram and redistribute the excess pixels in each bin
avgBinIncr = floor(totalExcess/numBins);
upperLimit = clipLimit - avgBinIncr; % bins larger than this will be
                                     % set to clipLimit

% this loop should speed up the operation by putting multiple pixels
% into the "obvious" places first
for k=1:numBins
  if imgHist(k) > clipLimit
    imgHist(k) = clipLimit;
  else
    if imgHist(k) > upperLimit % high bin count
      totalExcess = totalExcess - (clipLimit - imgHist(k));
      imgHist(k) = clipLimit;
    else
      totalExcess = totalExcess - avgBinIncr;
      imgHist(k) = imgHist(k) + avgBinIncr;      
    end
  end
end
 
% this loops redistributes the remaining pixels, one pixel at a time
k = 1;
while (totalExcess ~= 0)
  %keep increasing the step as fewer and fewer pixels remain for
  %the redistribution (spread them evenly)
  stepSize = max(floor(numBins/totalExcess),1);
  for m=k:stepSize:numBins
    if imgHist(m) < clipLimit
      imgHist(m) = imgHist(m)+1;
      totalExcess = totalExcess - 1; %reduce excess
      if totalExcess == 0
        break;
      end
    end
  end
  
  k = k+1; %prevent from always placing the pixels in bin #1
  if k > numBins % start over if numBins was reached
    k = 1;
  end
end

CLAHE和AHE中另一個重要的問題：插值。

將圖像進行分塊處理，若每塊中的像素點僅通過該塊中的映射函數進行變換，則會導致最終圖像呈塊狀效應：

為了解決這個問題，我們需要利用插值運算，也就是每個像素點出的值由它周圍4個子塊的映射函數值進行雙線性插值得到，如下圖：

上圖中，為了求藍色像素點處的值，需要利用它周圍四個子塊的映射函數分別做變換得到四個映射值，再對這四個值做雙線性插值即可。

當 然對於邊界處的像素點則不是通過四個子塊進行插值，如上圖紅色像素點直接以一個子塊的映射函數做變換，綠色像素則以兩個子塊做映射函數做線性插值。這里講 的邊界處像素是指落在圖像左上角，左下角、右上角，右下角的四個子塊中心像素點圍成的四邊形之外的像素。如下圖，將圖像分為8x8子塊，邊界像素即落在灰 色區域的像素點。

利用插值對圖像進行處理的整體架構如下：

 

for (Y = 0; Y <= TileY; Y++) //TileY為Y方向網格數
        {
            if (Y == 0) 
            {  
                SubY = TileYDim >> 1;  YU = 0; YB = 0;
            }
            else if (Y == TileY) 
            {   
                SubY = TileYDim >> 1;    YU = TileY-1;  YB = YU;
            }
            else 
            {  
                SubY = TileYDim; YU = Y - 1; YB = Y;
            }
            for (X = 0; X <= TileX; X++) //TileX為X方向網格數
            {
                if (X == 0) 
                {   
                    SubX = TileXDim >> 1; XL = 0; XR = 0;
                }
                else if (X == TileX)
                {   
                    SubX = TileXDim >> 1;  XL = TileX - 1; XR = XL;
                }
                else 
                {   
                    SubX = TileXDim; XL = X - 1; XR = X;
                }
                MapLU = &pMapArray[numBins * (YU * TileX + XL)];//左上角映射函數
                MapRU = &pMapArray[numBins * (YU * TileX + XR)];//右上角映射函數
                MapLB = &pMapArray[numBins * (YB * TileX + XL)];//左下角映射函數
                MapRB = &pMapArray[numBins * (YB * TileX + XR)];//右下角映射函數        
                Interpolate(pImPointer,Stride,Channel,MapLU,MapRU,MapLB,MapRB,SubX,SubY,aLUT);//插值
                pImPointer += SubX;          
            }
            pImPointer+=(SubY-1)*Stride;
        }

 

注意的是，上述循環需要（TileX+1）*（TileY+1）次，而不是TileX*TileY次。，原因很簡單，以X方向為例，兩側邊界處的半寬子塊（灰色區）也各需要處理一次，如下圖：

通過雙線性插值可以基本消除塊狀效應：

對 於彩色圖像，三通處理分開處理會導致嚴重的偏色，故我們可以將其進行顏色空間轉換（如RGB轉為HSV），然后僅對亮度分量處理，再反變換回RGB空間。 不過網上有高手將R、G、B統一處理[2]（也就相當於把一個像素點拆成三個像素點），這樣得到的效果也不錯，而且省去了顏色空間轉換的時間，我們這里也 仿照他來吧：

另外，CLAHE對霧天圖像處理效果也不錯：

至於編程，我基本就是翻譯adaphisteq.m，另外還有一些參考資源：CLAHE代碼（MATLAB）

這里給出編譯好的文件，有興趣的朋友可以下載看看：ImageProcess(CLAHE)

參考：

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

[2]http://www.cnblogs.com/Imageshop/archive/2013/04/07/3006334.html