自适应直方图均衡（AHE）和限制对比度的自适应直方图均衡（CLAHE）

本文翻译自https://en.wikipedia.org/wiki/Adaptive_histogram_equalization，如有错误还望海涵。。

自适应的直方图均衡（Adaptive Histogram Equalization）

1.算法简介

AHE是一种用来改善图像对比度的图像处理技术，它与传统的（普通）直方图均衡相比，不同点主要在于，AHE通过计算图像每一个显著区域的直方图，来重新分布图像的亮度值，因此它更适合于用来改善图像的局部对比度，以及增强图像边缘信息，利于分割。

但是，AHE有一个缺陷，就是他在增强对比度的同时也会增强图像同质（均匀）区域的噪声，因此作为AHE的改进，CLAHE可以有效降低这种噪声的增强。

2.算法的意义

传统的直方图均衡在对整幅图像进行变换时都采用相同的直方图变换，这种方法对于那些像素分布均匀的图像来说效果很好，但是对于那种包含明显较亮或较暗区域的图像来说，往往不能起到显著的增强效果。

AHE通过计算每一个像素邻域的变换函数来对每个像素执行直方图增强，该算法最早被用在航天器驾驶舱的显示图像增强上。其最简单的形式是，基于该像素的方形邻域的直方图来均衡化每一个像素，如下图，这种基于直方图进行增强的思想其实跟普通的直方图增强完全相同，因为这种变换函数与像素邻域的累积分布函数（CDF）是相称的。

对于图像的边界像素应该特殊处理，因为这些像素的邻域并不完全包含在图像内，比如上图中蓝色框，这可以通过像素行或列的镜像来解决；直接对边界像素行进行复制是不可取的，因为这样会造成像素邻域直方图出现很高的峰值。

3.AHE算法的性质

（1）邻域的size是一个邻域长度的尺度参数，当这个参数较大时，会降低对比度；反之增强对比度；

（ 2）由于直方图增强的性质，AHE中的像素变换结果是跟其在邻域像素中的顺序排列成正比的，因此我们可以利用专门的硬件来实现这种邻域中心像素和其他像素之间的比较；一种非归一化的计算方法，可以通过对每一个像素加上一个比中心像素小的值来获得，或者在每一个像素上加上一个邻域像素均值。

（3）如果一个像素邻域是同质的或均匀的，则该邻域的直方图可能会尖状化，即出现很高的峰值，这时变换函数就会把一个很窄的像素值范围映射到变换结果的整个像素值范围，这也正是上面提到的AHE为什么会增强均匀像素邻域噪声的原因所在。

限制对比度的自适应直方图增强（Contrast Limited AHE）

1.基本原理

CLAHE与AHE最大的不同在于前者对对比度进行了限制，这一特性也可以被应用到全局的直方图均衡中，即Contrast Limited HE，简称CLHE，但实际中它很少被用到。CLAHE中，每一个像素邻域都要进行对比度限制，从而得到对应的变换函数，被用来降低AHE中噪声的增强，这主要是通过限制AHE中的对比度增强来实现的。像素周围邻域噪声的增强主要是由变换函数的斜率造成的，由于像素邻域的噪声与邻域的CDF成正比，因此也与邻域直方图在该中心像素位置的值成正比，CLAHE之所以能够限制对比度，是因为它在计算邻域的CDF之前在指定阈值处对直方图进行了修剪，如下图所示，这一做法不仅限制了CDF的斜率，也限制了变换函数的斜率，其中对直方图进行切割所使用的阈值，被称作修剪限制度（clip
limit），这个参数不仅依赖于直方图的归一化，而且依赖于像素邻域的size大小，通常设为3到4之间。

上图中，并没有将修剪掉的那部分直方图直接扔掉，而是采用了更妥当的办法，就是将这些被修剪掉的部分重新均匀的分布到直方图中，从而生成新的直方图。

另外，在这种重新分布下，有可能导致部分区域会再次超过clip limit，如上图右边的图中标记为绿色的那部分区域，这时我们需要调大clip limit以达到实际需要，如果仍然不能让你满意，可以尝试使用递归的过程来重复这一重新分布直方图的步骤，直到超出的部分对我们的结果影响可以忽略不计。

2.插值加速

上面提到的自适应的直方图均衡算法，不论有没有对对比度进行限制，都需要对每一个像素的像素邻域计算邻域直方图和对应的变换函数，这个过程是相当耗时的。

因此，在保证算法效果的基础上，插值算法被用来提升计算效率，具体做法如下：

-将图像分为相同size的矩形tile，如下图右，通常分为8x8，即64个tiles；

-分别计算每一个tile的直方图、CDF以及变换函数；对于tiles的中心像素（下图左侧的黑色矩形块），由于满足于上述求得的变换函数，因此可以直接使用所在tile的变换函数得到，而其他像素则需要通过与其最邻近的至少四个像素所在的tiles（即四邻域，还可以8邻域）的变换函数进行插值得到。其中，下图左中蓝色区域标记的像素点是用双线性插值得到（注，这部分像素点占一幅图像的大部分）；边界部分的像素点（绿色标记的区域）是用线性插值得到的；而靠近图像角点部分（红色标记）的像素点则是通过角点所在tile的变换函数得到。

其中插值系数反映了距离最近的tiles的中心像素之间的位置，因此这个插值过程是连续的，因为插值得到的像素点可以近似作中心像素点。

这一过程虽然增加了线性插值的计算量，但是却大大降低了变换函数的计算次数，所以提高了计算速度。

3.CLAHE的源码参考

/*
 * ANSI C code from the article
 * "Contrast Limited Adaptive Histogram Equalization"
 * by Karel Zuiderveld, [email protected]
 * in "Graphics Gems IV", Academic Press, 1994
 *
 *
 *  These functions implement Contrast Limited Adaptive Histogram Equalization.
 *  The main routine (CLAHE) expects an input image that is stored contiguously in
 *  memory;  the CLAHE output image overwrites the original input image and has the
 *  same minimum and maximum values (which must be provided by the user).
 *  This implementation assumes that the X- and Y image resolutions are an integer
 *  multiple of the X- and Y sizes of the contextual regions. A check on various other
 *  error conditions is performed.
 *
 *  #define the symbol BYTE_IMAGE to make this implementation suitable for
 *  8-bit images. The maximum number of contextual regions can be redefined
 *  by changing uiMAX_REG_X and/or uiMAX_REG_Y; the use of more than 256
 *  contextual regions is not recommended.
 *
 *  The code is ANSI-C and is also C++ compliant.
 *
 *  Author: Karel Zuiderveld, Computer Vision Research Group,
 *         Utrecht, The Netherlands ([email protected])
 */

#ifdef BYTE_IMAGE
typedef unsigned char kz_pixel_t;     /* for 8 bit-per-pixel images */
#define uiNR_OF_GREY (256)
#else
typedef unsigned short kz_pixel_t;     /* for 12 bit-per-pixel images (default) */
# define uiNR_OF_GREY (4096)
#endif

/******** Prototype of CLAHE function. Put this in a separate include file. *****/
int CLAHE(kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes, kz_pixel_t Min,
      kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,
      unsigned int uiNrBins, float fCliplimit);

/*********************** Local prototypes ************************/
static void ClipHistogram (unsigned long*, unsigned int, unsigned long);
static void MakeHistogram (kz_pixel_t*, unsigned int, unsigned int, unsigned int,
        unsigned long*, unsigned int, kz_pixel_t*);
static void MapHistogram (unsigned long*, kz_pixel_t, kz_pixel_t,
           unsigned int, unsigned long);
static void MakeLut (kz_pixel_t*, kz_pixel_t, kz_pixel_t, unsigned int);
static void Interpolate (kz_pixel_t*, int, unsigned long*, unsigned long*,
    unsigned long*, unsigned long*, unsigned int, unsigned int, kz_pixel_t*);

/**************     Start of actual code **************/
#include <stdlib.h>             /* To get prototypes of malloc() and free() */

const unsigned int uiMAX_REG_X = 16;      /* max. # contextual regions in x-direction */
const unsigned int uiMAX_REG_Y = 16;      /* max. # contextual regions in y-direction */

/************************** main function CLAHE ******************/
int CLAHE (kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes,
     kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,
          unsigned int uiNrBins, float fCliplimit)
/*   pImage - Pointer to the input/output image
 *   uiXRes - Image resolution in the X direction
 *   uiYRes - Image resolution in the Y direction
 *   Min - Minimum greyvalue of input image (also becomes minimum of output image)
 *   Max - Maximum greyvalue of input image (also becomes maximum of output image)
 *   uiNrX - Number of contextial regions in the X direction (min 2, max uiMAX_REG_X)
 *   uiNrY - Number of contextial regions in the Y direction (min 2, max uiMAX_REG_Y)
 *   uiNrBins - Number of greybins for histogram ("dynamic range")
 *   float fCliplimit - Normalized cliplimit (higher values give more contrast)
 * The number of "effective" greylevels in the output image is set by uiNrBins; selecting
 * a small value (eg. 128) speeds up processing and still produce an output image of
 * good quality. The output image will have the same minimum and maximum value as the input
 * image. A clip limit smaller than 1 results in standard (non-contrast limited) AHE.
 */
{
    unsigned int uiX, uiY;          /* counters */
    unsigned int uiXSize, uiYSize, uiSubX, uiSubY; /* size of context. reg. and subimages */
    unsigned int uiXL, uiXR, uiYU, uiYB;  /* auxiliary variables interpolation routine */
    unsigned long ulClipLimit, ulNrPixels;/* clip limit and region pixel count */
    kz_pixel_t* pImPointer;           /* pointer to image */
    kz_pixel_t aLUT[uiNR_OF_GREY];        /* lookup table used for scaling of input image */
    unsigned long* pulHist, *pulMapArray; /* pointer to histogram and mappings*/
    unsigned long* pulLU, *pulLB, *pulRU, *pulRB; /* auxiliary pointers interpolation */

    if (uiNrX > uiMAX_REG_X) return -1;       /* # of regions x-direction too large */
    if (uiNrY > uiMAX_REG_Y) return -2;       /* # of regions y-direction too large */
    if (uiXRes % uiNrX) return -3;      /* x-resolution no multiple of uiNrX */
    if (uiYRes & uiNrY) return -4;      /* y-resolution no multiple of uiNrY */
    if (Max >= uiNR_OF_GREY) return -5;       /* maximum too large */
    if (Min >= Max) return -6;          /* minimum equal or larger than maximum */
    if (uiNrX < 2 || uiNrY < 2) return -7;/* at least 4 contextual regions required */
    if (fCliplimit == 1.0) return 0;      /* is OK, immediately returns original image. */
    if (uiNrBins == 0) uiNrBins = 128;      /* default value when not specified */

    pulMapArray=(unsigned long *)malloc(sizeof(unsigned long)*uiNrX*uiNrY*uiNrBins);
    if (pulMapArray == 0) return -8;      /* Not enough memory! (try reducing uiNrBins) */

    uiXSize = uiXRes/uiNrX; uiYSize = uiYRes/uiNrY;  /* Actual size of contextual regions */
    ulNrPixels = (unsigned long)uiXSize * (unsigned long)uiYSize;

    if(fCliplimit > 0.0) {          /* Calculate actual cliplimit     */
       ulClipLimit = (unsigned long) (fCliplimit * (uiXSize * uiYSize) / uiNrBins);
       ulClipLimit = (ulClipLimit < 1UL) ? 1UL : ulClipLimit;
    }
    else ulClipLimit = 1UL<<14;          /* Large value, do not clip (AHE) */
    MakeLut(aLUT, Min, Max, uiNrBins);      /* Make lookup table for mapping of greyvalues */
    /* Calculate greylevel mappings for each contextual region */
    for (uiY = 0, pImPointer = pImage; uiY < uiNrY; uiY++) {
    for (uiX = 0; uiX < uiNrX; uiX++, pImPointer += uiXSize) {
        pulHist = &pulMapArray[uiNrBins * (uiY * uiNrX + uiX)];
        MakeHistogram(pImPointer,uiXRes,uiXSize,uiYSize,pulHist,uiNrBins,aLUT);
        ClipHistogram(pulHist, uiNrBins, ulClipLimit);
        MapHistogram(pulHist, Min, Max, uiNrBins, ulNrPixels);
    }
    pImPointer += (uiYSize - 1) * uiXRes;          /* skip lines, set pointer */
    }

    /* Interpolate greylevel mappings to get CLAHE image */
    for (pImPointer = pImage, uiY = 0; uiY <= uiNrY; uiY++) {
    if (uiY == 0) {                      /* special case: top row */
        uiSubY = uiYSize >> 1;  uiYU = 0; uiYB = 0;
    }
    else {
        if (uiY == uiNrY) {                  /* special case: bottom row */
        uiSubY = uiYSize >> 1;    uiYU = uiNrY-1;     uiYB = uiYU;
        }
        else {                      /* default values */
        uiSubY = uiYSize; uiYU = uiY - 1; uiYB = uiYU + 1;
        }
    }
    for (uiX = 0; uiX <= uiNrX; uiX++) {
        if (uiX == 0) {                  /* special case: left column */
        uiSubX = uiXSize >> 1; uiXL = 0; uiXR = 0;
        }
        else {
        if (uiX == uiNrX) {              /* special case: right column */
            uiSubX = uiXSize >> 1;  uiXL = uiNrX - 1; uiXR = uiXL;
        }
        else {                      /* default values */
            uiSubX = uiXSize; uiXL = uiX - 1; uiXR = uiXL + 1;
        }
        }

        pulLU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXL)];
        pulRU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXR)];
        pulLB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXL)];
        pulRB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXR)];
        Interpolate(pImPointer,uiXRes,pulLU,pulRU,pulLB,pulRB,uiSubX,uiSubY,aLUT);
        pImPointer += uiSubX;              /* set pointer on next matrix */
    }
    pImPointer += (uiSubY - 1) * uiXRes;
    }
    free(pulMapArray);                      /* free space for histograms */
    return 0;                          /* return status OK */
}
void ClipHistogram (unsigned long* pulHistogram, unsigned int
             uiNrGreylevels, unsigned long ulClipLimit)
/* This function performs clipping of the histogram and redistribution of bins.
 * The histogram is clipped and the number of excess pixels is counted. Afterwards
 * the excess pixels are equally redistributed across the whole histogram (providing
 * the bin count is smaller than the cliplimit).
 */
{
    unsigned long* pulBinPointer, *pulEndPointer, *pulHisto;
    unsigned long ulNrExcess, ulUpper, ulBinIncr, ulStepSize, i;
    long lBinExcess;

    ulNrExcess = 0;  pulBinPointer = pulHistogram;
    for (i = 0; i < uiNrGreylevels; i++) { /* calculate total number of excess pixels */
    lBinExcess = (long) pulBinPointer[i] - (long) ulClipLimit;
    if (lBinExcess > 0) ulNrExcess += lBinExcess;      /* excess in current bin */
    };

    /* Second part: clip histogram and redistribute excess pixels in each bin */
    ulBinIncr = ulNrExcess / uiNrGreylevels;          /* average binincrement */
    ulUpper =  ulClipLimit - ulBinIncr;     /* Bins larger than ulUpper set to cliplimit */

    for (i = 0; i < uiNrGreylevels; i++) {
      if (pulHistogram[i] > ulClipLimit) pulHistogram[i] = ulClipLimit; /* clip bin */
      else {
      if (pulHistogram[i] > ulUpper) {        /* high bin count */
          ulNrExcess -= pulHistogram[i] - ulUpper; pulHistogram[i]=ulClipLimit;
      }
      else {                    /* low bin count */
          ulNrExcess -= ulBinIncr; pulHistogram[i] += ulBinIncr;
      }
       }
    }

    while (ulNrExcess) {   /* Redistribute remaining excess  */
    pulEndPointer = &pulHistogram[uiNrGreylevels]; pulHisto = pulHistogram;

    while (ulNrExcess && pulHisto < pulEndPointer) {
        ulStepSize = uiNrGreylevels / ulNrExcess;
        if (ulStepSize < 1) ulStepSize = 1;          /* stepsize at least 1 */
        for (pulBinPointer=pulHisto; pulBinPointer < pulEndPointer && ulNrExcess;
         pulBinPointer += ulStepSize) {
        if (*pulBinPointer < ulClipLimit) {
            (*pulBinPointer)++;     ulNrExcess--;      /* reduce excess */
        }
        }
        pulHisto++;          /* restart redistributing on other bin location */
    }
    }
}
void MakeHistogram (kz_pixel_t* pImage, unsigned int uiXRes,
        unsigned int uiSizeX, unsigned int uiSizeY,
        unsigned long* pulHistogram,
        unsigned int uiNrGreylevels, kz_pixel_t* pLookupTable)
/* This function classifies the greylevels present in the array image into
 * a greylevel histogram. The pLookupTable specifies the relationship
 * between the greyvalue of the pixel (typically between 0 and 4095) and
 * the corresponding bin in the histogram (usually containing only 128 bins).
 */
{
    kz_pixel_t* pImagePointer;
    unsigned int i;

    for (i = 0; i < uiNrGreylevels; i++) pulHistogram[i] = 0L; /* clear histogram */

    for (i = 0; i < uiSizeY; i++) {
    pImagePointer = &pImage[uiSizeX];
    while (pImage < pImagePointer) pulHistogram[pLookupTable[*pImage++]]++;
    pImagePointer += uiXRes;
    pImage = &pImagePointer[-uiSizeX];
    }
}

void MapHistogram (unsigned long* pulHistogram, kz_pixel_t Min, kz_pixel_t Max,
           unsigned int uiNrGreylevels, unsigned long ulNrOfPixels)
/* This function calculates the equalized lookup table (mapping) by
 * cumulating the input histogram. Note: lookup table is rescaled in range [Min..Max].
 */
{
    unsigned int i;  unsigned long ulSum = 0;
    const float fScale = ((float)(Max - Min)) / ulNrOfPixels;
    const unsigned long ulMin = (unsigned long) Min;

    for (i = 0; i < uiNrGreylevels; i++) {
    ulSum += pulHistogram[i]; pulHistogram[i]=(unsigned long)(ulMin+ulSum*fScale);
    if (pulHistogram[i] > Max) pulHistogram[i] = Max;
    }
}

void MakeLut (kz_pixel_t * pLUT, kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrBins)
/* To speed up histogram clipping, the input image [Min,Max] is scaled down to
 * [0,uiNrBins-1]. This function calculates the LUT.
 */
{
    int i;
    const kz_pixel_t BinSize = (kz_pixel_t) (1 + (Max - Min) / uiNrBins);

    for (i = Min; i <= Max; i++)  pLUT[i] = (i - Min) / BinSize;
}

void Interpolate (kz_pixel_t * pImage, int uiXRes, unsigned long * pulMapLU,
     unsigned long * pulMapRU, unsigned long * pulMapLB,  unsigned long * pulMapRB,
     unsigned int uiXSize, unsigned int uiYSize, kz_pixel_t * pLUT)
/* pImage      - pointer to input/output image
 * uiXRes      - resolution of image in x-direction
 * pulMap*     - mappings of greylevels from histograms
 * uiXSize     - uiXSize of image submatrix
 * uiYSize     - uiYSize of image submatrix
 * pLUT           - lookup table containing mapping greyvalues to bins
 * This function calculates the new greylevel assignments of pixels within a submatrix
 * of the image with size uiXSize and uiYSize. This is done by a bilinear interpolation
 * between four different mappings in order to eliminate boundary artifacts.
 * It uses a division; since division is often an expensive operation, I added code to
 * perform a logical shift instead when feasible.
 */
{
    const unsigned int uiIncr = uiXRes-uiXSize; /* Pointer increment after processing row */
    kz_pixel_t GreyValue; unsigned int uiNum = uiXSize*uiYSize; /* Normalization factor */

    unsigned int uiXCoef, uiYCoef, uiXInvCoef, uiYInvCoef, uiShift = 0;

    if (uiNum & (uiNum - 1))   /* If uiNum is not a power of two, use division */
    for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
     uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {
    for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
         uiXCoef++, uiXInvCoef--) {
        GreyValue = pLUT[*pImage];           /* get histogram bin value */
        *pImage++ = (kz_pixel_t ) ((uiYInvCoef * (uiXInvCoef*pulMapLU[GreyValue]
                      + uiXCoef * pulMapRU[GreyValue])
                + uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
                      + uiXCoef * pulMapRB[GreyValue])) / uiNum);
    }
    }
    else {               /* avoid the division and use a right shift instead */
    while (uiNum >>= 1) uiShift++;           /* Calculate 2log of uiNum */
    for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
         uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {
         for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
           uiXCoef++, uiXInvCoef--) {
           GreyValue = pLUT[*pImage];      /* get histogram bin value */
           *pImage++ = (kz_pixel_t)((uiYInvCoef* (uiXInvCoef * pulMapLU[GreyValue]
                      + uiXCoef * pulMapRU[GreyValue])
                + uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
                      + uiXCoef * pulMapRB[GreyValue])) >> uiShift);
        }
    }
    }
}