对比度受限的自适应直方图均衡化(CLAHE)

直方图均衡化(HE)简介

图像的灰度级

  在计算机中,一幅图像由若干个像素组成,这些像素有自己的位置和像素值.

  一幅彩色图由RGB三个通道的像素组成,一副灰度图由单通道的像素组成.

  在一副灰度图中,像素值的取值范围即为灰度级,在计算机中,灰度级为256,即[0, 255],灰度级为离散值.

图像的灰度直方图

  图像的直方图就是对图像的每个灰度级中有多少个像素点的统计,如下图:

直方图均衡化(HE)

  显然,对于一些图像而言,它的像素值可能集中分布在某一块区域,那么像素点与像素点之间的色差就会很小,这就会导致像素点之间的对比度很小,图像中的一些细节就会不清晰.

  通过直方图均衡化对图像进行处理后,使原始图像的像素值尽可能均匀的分布在所有灰度级中(0~255),这样就拉抻了图像的色域,像素点之间的对比度得到了增强,图像的细节也就更为明显.

HE的数学原理

  显然,需要找到一组映射关系T,将原图像的像素值x映射到目标图像中对应的位置,且值为T(x),使目标图像的直方图均匀的分布在整个灰度级中.

  根据任务目标可以分析出,是为了能让目标图像的像素点均匀分布在整个灰度区间而去寻找的映射关系T,因此要从原始图像和目标图像的直方图入手,设:
  图像的灰度级为L;
  原始图像的像素值为x,则x∈[0,L-1],假设x为连续值;
  像素数量为N;
  像素值为x的像素点数量为Nx;

则原始图像的直方图为hist(x)=Nx.
设像素值的频率为概率,则原始图像的概率分布为F(x)=Nx/N,且
  设目标图像的概率分布为H(s),h(s)为目标图像的概率密度函数,f(x)为原始图像的概率密度函数,s=T(x),则H(s)=F(T(x)).
等式两边都对s求导得:

  因为目标图像的像素值应均匀分布在[0,L-1]中,所以目标图像的概率密度函数h(s)=1,所以ds/dx=f(x),即

对等式两边分别求积分得T(x)=F(x).
  显然,为了使目标图像的像素点均匀分布在整个灰度区间而去寻找的映射关系T(x)就是原始图像的概率分布函数F(x),然而像素值x实际上为离散值,因此T(x)为原始图像的累计分布函数.

自适应直方图均衡化(AHE)

  然而,HE是对图像全局进行处理的,因此对于一些原本颜色偏黑的图像,将色域拉抻后得到的图像会发白,并且局部的对比度不能有效提高.

  为了提高图像的局部对比度,有人提出了AHE,即将图像划分为若干子块,对子块进行HE处理.

  对于这种处理方式,不难想到,AHE对局部对比度提高过大,会导致图像失真.此外,局部对比度过高还会放大图像中的噪声.

对比度受限的自适应直方图均衡化(CLAHE)

  对子块中得到的直方图进行剪裁,使其幅值低于给定的上限.由于剪裁掉的部分是像素点的数量,因此不能丢掉,所以要将这部分像素点均匀的分布在整个灰度区间上,以保证没有像素点缺失,如下图:

  可以看到,这时直方图又会整体上升了一个高度,会超过我们设置的上限.其实在具体实现的时候有很多解决方法,你可以多重复几次裁剪过程,使得上升的部分变得微不足道,或是用另一种常用的方法:

  设裁剪值为ClipLimit,求直方图中高于该值的部分的和totalExcess,此时假设将totalExcess均分给所有灰度级,求出这样导致的直方图整体上升的高度L=totalExcess/N,以upper= ClipLimit-L为界限对直方图进行如下处理：

  (1)若幅值高于ClipLimit，直接置为ClipLimit；
  (2)若幅值处于Upper和ClipLimit之间，将其填补至ClipLimit；
  (3)若幅值低于Upper，直接填补L个像素点；

经过上述操作,用来填补的像素点个数通常会略小于totalExcess,也就是还有一些剩余的像素点没分出去.这时我们可以再把这些点均匀地分给那些目前幅值仍然小于ClipLimit的灰度值.

代码如下:

function imgHist = clipHistogram(imgHist, clipLimit, numBins)
%求出原图像子块的对比度受限的直方图
totalExcess = sum(max(imgHist - clipLimit,0)); 
avgBinIncr = floor(totalExcess/numBins);
upperLimit = clipLimit - avgBinIncr;

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

k = 1;
while (totalExcess ~= 0)
  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
  if stepSize~=1
    k = k+1; 
    if k > numBins 
      k = 1;
    end
  end
end
end

得到每一块的映射关系T(x)

function mapping = makeMapping(imgHist,  fullRange,  numPixInTile)

histSum = cumsum(imgHist);
valSpread  = fullRange(2) - fullRange(1);
  scale =  valSpread/numPixInTile;
  mapping = min(fullRange(1) + histSum*scale, fullRange(2)); 
end

AHE与CLAHE中存在的另一个问题

  将图像进行分块处理，若每块中的像素点仅通过该块中的映射函数进行变换，则会导致最终图像呈块状,如下图：

代码如下:

function claheI = makeClaheImage(I, tileMappings, numTiles, dimTile)
% 未进行插值
claheI = I;
claheI(:) = 0;
imgTileRow=1;
for i = 1:numTiles(1)
    imgTileCol=1;
    for j = 1:numTiles(2)
        for row = 1:dimTile(1)     
            for col = 1: dimTile(2)
                rowidx = imgTileRow+row-1;
                colidx = imgTileCol+col-1;
                claheI(rowidx, colidx) = ...
                    tileMappings{i, j}(I(rowidx,colidx)+1);
            end
        end
        imgTileCol = imgTileCol + dimTile(1);
    end
    imgTileRow = imgTileRow + dimTile(2);
end
end

插值

  为了解决最终图像呈块状的效应,采用插值的方式计算最终的像素值.
  将原图像分块,以8*8为例,如下图:

  划分边界,如下图:

  遍历顺序如下图:

  显然,遍历顺序不是划分的8*8的块,而是错开遍历,这样每次遍历的块中就包含了原始划分规则中的边界,进行插值即可消除块状效应.

大体规则为:

  (1)边角位:不插值,直接将像素值带入Ti(x);

  (2)最上与最下的一行:将像素值带入其所在子域与右邻子域的映射函数,然后做线性插值;

  (3)最左与最右的一列:将像素值带入其所在子域与上邻子域的映射函数,然后做线性插值;

  (4)常规位置:将像素值带入其所在子域与左,上,左上子域的映射函数,然后做双线性插值.

代码如下:

function claheI = makeClaheImage(I, tileMappings, numTiles, dimTile)

claheI = I;
claheI(:) = 0;

imgTileRow=1;
for k=1:numTiles(1)+1
  if k == 1  
    imgTileNumRows = dimTile(1)/2;
    mapTileRows = [1 1];
  else 
    if k == numTiles(1)+1  
      imgTileNumRows = dimTile(1)/2;
      mapTileRows = [numTiles(1) numTiles(1)];
    else 
      imgTileNumRows = dimTile(1); 
      mapTileRows = [k-1, k]; 
    end
  end
  
  imgTileCol=1;
  for l=1:numTiles(2)+1
    if l == 1 
      imgTileNumCols = dimTile(2)/2;
      mapTileCols = [1, 1];
    else
      if l == numTiles(2)+1 
        imgTileNumCols = dimTile(2)/2;
        mapTileCols = [numTiles(2), numTiles(2)];
      else 
        imgTileNumCols = dimTile(2);
        mapTileCols = [l-1, l]; 
      end
    end
    
    ulMapTile = tileMappings{mapTileRows(1), mapTileCols(1)};
    urMapTile = tileMappings{mapTileRows(1), mapTileCols(2)};
    blMapTile = tileMappings{mapTileRows(2), mapTileCols(1)};
    brMapTile = tileMappings{mapTileRows(2), mapTileCols(2)};
    
    normFactor = imgTileNumRows*imgTileNumCols; 

    for i = 1:imgTileNumRows
        
        for j = 1: imgTileNumCols
            rowidx = imgTileRow+i-1;
            colidx  = imgTileCol+j-1;
            
            claheI(rowidx, colidx) = ...
                ((imgTileNumRows-i)*((imgTileNumCols-j)*double(ulMapTile(I(rowidx,colidx)+1))+...
                    j*double(urMapTile(I(rowidx,colidx)+1)))+...
                  i*((imgTileNumCols-j)*double(blMapTile(I(rowidx,colidx)+1))+...
                    j*double(brMapTile(I(rowidx,colidx)+1))))...
                 /normFactor;                
        end
    end
    imgTileCol = imgTileCol + imgTileNumCols;    
  end
  imgTileRow = imgTileRow + imgTileNumRows;
end 
end

原图像与最终结果对比

原始图像:

最终结果:

原始图像直方图:

最终图像直方图:

完整代码:

  https://github.com/Qyokizzzz/AI-Algorithm/tree/master/CLAHE

参考:

  https://blog.csdn.net/superjunenaruto/article/details/52431941
  https://www.cnblogs.com/luo-peng/p/4930013.html