CLAHE算法

直方图均衡（HE）

直方图均衡 histogram equalization是图像增强的一种方法。这个过程的变换函数为：
$$s=T(r)=(L-1){\int }_0^r{p}_r(w)\mathrm{d}w(\star )$$

这个公式的来源可以参考 Gonzalez 的数字图像处理教材，这里简单推导一下：

• 如果我们采用$(\star )$式的变换函数，则有：
  $$\frac{\mathrm{d}s}{\mathrm{d}r}=\frac{\mathrm{d}T(r)}{\mathrm{d}r}=(L-1)\frac{\mathrm{d}}{\mathrm{d}r}\left[{\int }_0^r{p}_r(w)\mathrm{d}w\right]=(L-1)p_r(r)$$

• 而同时在将灰度$r$映射到灰度$s$过程中对于$\frac{\mathrm{d}s}{\mathrm{d}r}$又有：

$$\begin{gathered} {\int }_0^s{P}_s(w)\mathrm{d}w={\int }_0^r{P}_r(w)\mathrm{d}w \\ \Rightarrow p_s(s)=\frac{\mathrm{d}({\int }_0^r{P}_r(w)\mathrm{d}w)}{\mathrm{d}s}=p_r(r)\frac{\mathrm{d}r}{\mathrm{d}s} \end{gathered}$$

• 结合上述两个式子可以得到：

$$p_s(s)=p_r(r)\frac{\mathrm{d}r}{\mathrm{d}s}=p_r(r)\frac{1}{(L-1)p_r(r)}=\frac{1}{L-1},0⩽s⩽L-1$$

• 也即一个不均匀的灰度分布变成了一个完全均匀的灰度分布

当然上述的结论是对于一个完全连续分布的灰度直方图来说，实际中灰度是离散分布的，则公式变为：
$$s_k=T\left(r_k\right)=(L-1)\sum _{j=0}^k{p}_r\left(r_j\right)=\frac{(L-1)}{MN}\sum _{j=0}^k{n}_j(k=0,1,\cdots ,L-1)$$

算法改进（AHE/CLAHE）

貌似这个直方图均衡效果很好，但实际不然，从它的两个缺点衍生出了AHE和CLAHE算法：

参考文献：Zuiderveld K. Contrast limited adaptive histogram equalization[J]. Graphics gems, 1994: 474-485.

局部处理（双线性插值）

传统的直方图均衡是对于整个图像的，但实际上对整张图像直接进行直方图均衡是很不合理的，在图像中各个部分有一定差异的情况下进行全局的直方图均衡效果不好。所以产生了自适应直方图均衡算法（AHE），可以想到的是：

• 对于像素$A$（灰度值为$r_A$），对其本身及其邻域进行直方图均衡并获取映射函数$s=T_A(r)$，利用这个映射对像素A进行处理$s_A=T_A(r_A)$；

• 由于上述方法存在计算量太大的问题，我们可以采用另一种方法，也即对图像进行分块，对每一块分别进行直方图均衡；

• 对图像分块并对每一块进行直方图均衡看似降低了运算量，但可以预见其大概率会导致块效应的出现，故我们需要采用双线性插值来在一定程度上降低块效应。（双线性插值就是线性插值的推广，具体百度即可）。

  如图所示，对于白色像素点$O$点，灰度值为$r_O$，其周围最近的四个块的的灰度映射函数如图所示，则双线性插值的结果就是：
  $$s_O=(1-x)(1-y)g_A(r_O)+x(1-y)g_B(r_O)+(1-x)y{g}_C(r_O)+xy{g}_D(r_O)$$
  

对比度限制

直方图均衡还有一个缺点。举个极端的例子，一个完全黑色的图像在进行直方图均衡后会变成一个完全白色的图像，这是灰度映射中映射函数完全由灰度值的频数决定引起的。所以我们需要削弱高频数的灰度值的频数，以减轻这种不良效果。

即：获取图像直方图$\to$对直方图进行对比度限制$\to$对对比度限制后的直方图进行直方图均衡获取映射函数$\to$利用映射函数处理原图像

具体阈值处理及切分过程可以参考matlab源码（这是我基于源码编写的c++代码）：

基于直方图的图像增强算法 - 推杯问盏 - 博客园 (cnblogs.com)

int* ContrastLimit(int* hist)
{
    int* limithist = (int*)malloc(DEPTH * sizeof(int));
    int sum = 0;
    for (int i = 0; i < DEPTH; i++) {
        limithist[i] = hist[i];
        sum = sum + hist[i];
    }
    int clipLimit = CLIP_LIMIT_RATE * sum;
    //计算超出的像素的总个数及灰度值的个数
    int totalExcess = 0;
    int totalnum = 0;
    for (int i = 0; i < DEPTH; i++) {
        if (hist[i] > clipLimit) {
            totalExcess = totalExcess + hist[i] - clipLimit;
        }
    }
    //计算平均增加值（注意：由这部产生了多余值）
    int avgBinIncr = floor(totalExcess / (DEPTH - totalnum));
    //计算第二级阈值
    int upperLimit = clipLimit - avgBinIncr;
    //从灰度级1开始，对于每个灰度值进行如下判断：
    for (int i = 0; i < DEPTH; i++) {
        if (hist[i] > clipLimit) {
            limithist[i] = clipLimit;
        }
        else if (hist[i] > upperLimit) {
            //二级阈值直接设置为cliplimit；重新计算超出的像素个数
            totalExcess = totalExcess - (clipLimit - hist[i]);
            limithist[i] = clipLimit;
        }
        else {
            totalExcess = totalExcess - avgBinIncr;
            limithist[i] = hist[i] + avgBinIncr;
        }

    }
    //上述过程后仍会剩余有像素点未分配
    //设置步长来将剩余的像素值均分
    int step;
    while (totalExcess != 0) {
        step = max(DEPTH / totalExcess, 1);
        for (int k = 0;; k = k + step) {
            int i = k % 256;
            if (limithist[i] < clipLimit) {
                limithist[i] = limithist[i] + 1;
                totalExcess = totalExcess - 1;
            }
            if (totalExcess == 0) {
                break;
            }
        }
    }
    return limithist;
}

我的算法实现流程

由于需要自主完成所有过程，这是我理解的算法流程（处理灰度图像）：

• 将图像灰度化并调整大小，使其完全分割（调整大小可以采用镜像的方法补齐）

• 将图像分块，并计算每一块的直方图数组（二维数组）

• 对每一块的直方图进行对比度限制处理

• 利用限制后的直方图生成每一块的直方图均衡映射函数（二维数组）

• 判断像素点的位置并利用相应四个直方图映射函数进行双线性插值得到最终结果

  双线性插值需要注意的点：

  • 对于一个像素点，首先计算其所处块及其到中心的距离：

    x = i / GRID_SIZE;
    y = j / GRID_SIZE;
    xdiff = (double(i % GRID_SIZE) - center) / GRID_SIZE;
    ydiff = (double(j % GRID_SIZE) - center) / GRID_SIZE;
    

  • 根据像素坐标及差值的正负找到其临近的块（如图，红色部分有一个，绿色部分有两个，蓝色部分有四个，可以通过逻辑判断来实现），再根据距离获取权值。