计算机视觉算法探究：OpenCV CLAHE 算法详解| 社区征文

一、引言

2021 年 10 月开始学习 OpenCV 对比度受限的自适应直方图均衡 CLAHE，应用编程简单，了解详细算法却相当难。

创建 CLAHE 对象时，只传递了两个参数：clipLimit 和 tileGridSize，其中 clipLimit 是裁剪限制参数，tileGridSize 图像的分块个数。关于参数含义及相关的介绍请参考《OpenCV-Python自适应直方图均衡类CLAHE及方法详解》。

CLAHE 算法的基本步骤如下：

1. 将图像按参数 tileGridSize 切分为若干子块，这样图像就分成了 tileGridSize.height 行和 tileGridSize.width；
2. 对每个子块进行直方图均衡，计算直方图及累计直方图，得到每块原始灰度和均衡后的灰度映射表。关于直方图均衡请参考《数字图像处理：OpenCV直方图均衡算法研究及模拟实现》；
3. 对每个子块直方图的每个灰度级，使用 clipLimit 值做限定和裁剪；
4. 遍历输入图像每个点，以每点灰度值结合该点所在位置及周边关联分块进行灰度插值处理得到每点对应输出图像的灰度值。

看起来并不难，但在学习时查阅了各种公开资料，发现并不能解答学习时思考的一些问题，如：

1. 图像横向和纵向分块大小与图像的宽和高不能整除怎么处理？
2. CLIP 的剪裁是怎么实施的？
3. 插值处理具体算法怎样？

经过近 4 个多月断断续续性的学习，特别是对 OpenCV 自适应直方图均衡 CLAHE 源代码的深入解读，这些问题都得到了解决，下面就详细介绍一下。关于 OpenCV 自适应直方图均衡 CLAHE 详细源代码请参考《OpenCV自适应直方图均衡CLAHE C++源代码分享》。

二、图像分块坐标以及图像分块大小与图像的宽和高不能整除的处理

2.1、图像分块坐标

进行对比度受限的自适应直方图均衡处理时，首先是需要将图像按参数 tileGridSize 切分为若干子块，这样图像就分成了 tileGridSize.height 行和 tileGridSize.width 列。

对这种分块，每个分块在坐标体系进行标记的话，横向坐标范围为[0,tilesX_-1]，纵向坐标范围为[0,tilesY_]。这种标记分块的坐标，它与像素坐标存在映射关系，但是独立于像素坐标体系，老猿称这种分块的坐标为图像分块坐标。

2.2、不能整除的处理

当图像的宽（或高）不是对应横向（或纵向）分块数的整数倍时，老猿认为对于分块的处理有多种方式：

1. 将每个分块横向或纵向多加 1 个像素，最后一个分块的大小比前面分块小；
2. 将每个分块横向或纵向减去 1 个像素，最后一个分块的大小比前面分块大；
3. 将图像裁剪或补齐到可以整除的大小。

通过阅读源代码，OpenCV 中采用将图像补齐到可以整除的大小，即对于图像的宽（或高）不是对应横向（或纵向）分块数的整数倍时，将对应宽（或高）补齐到可以整除的最少像素素。

具体处理的源代码如下：

if (_src.size().width % tilesX_ == 0 && _src.size().height % tilesY_ == 0)
{
    tileSize = cv::Size(_src.size().width / tilesX_, _src.size().height / tilesY_);
    _srcForLut = _src;
}
else
{

    {
        cv::copyMakeBorder(_src, srcExt_, 0, tilesY_ - (_src.size().height % tilesY_), 0, 
                tilesX_ - (_src.size().width % tilesX_), cv::BORDER_REFLECT_101);
        tileSize = cv::Size(srcExt_.size().width / tilesX_, srcExt_.size().height / tilesY_);
        _srcForLut = srcExt_;
    }
}

上述代码中，src 是输入图像矩阵，tilesX_是横向分块数，tilesY_是纵向分块数，因此图像被分成了 tilesX_*tilesY_个分块。

三、CLIP 的赋值和裁剪过程

3.1、CLIP 的赋值过程

CLAHE 涉及 clipLimit 的关键源代码摘要如下：

    CLAHE_Impl::CLAHE_Impl(double clipLimit, int tilesX, int tilesY) :
        clipLimit_(clipLimit), tilesX_(tilesX), tilesY_(tilesY)
    {
    }

    void CLAHE_Impl::apply(cv::InputArray _src, cv::OutputArray _dst)
    {
	...
        int histSize = _src.type() == CV_8UC1 ? 256 : 65536;
        ...
         if (_src.size().width % tilesX_ == 0 && _src.size().height % tilesY_ == 0)
        {
            tileSize = cv::Size(_src.size().width / tilesX_, _src.size().height / tilesY_);
            _srcForLut = _src;
        }
        ...
        const int tileSizeTotal = tileSize.area();
        ...
        int clipLimit = 0;
        if (clipLimit_ > 0.0)
        {
            clipLimit = static_cast<int>(clipLimit_ * tileSizeTotal / histSize);
            clipLimit = std::max(clipLimit, 1);
        }
        ...
    }

以上代码就是 OpenCV 自适应直方图均衡 CLAHE 对应源代码中关于 clipLimit 赋值处理的相关代码，可以看到，类设置方法中对 clipLimit 设置后，其值会保存在类私有变量 clipLimit_ 中，最终进行 apply 自适应直方图均衡处理时，采用局部变量 clipLimit = clipLimit_ * tileSizeTotal / histSize，并取 clipLimit 和 1 中间的最大值。

可以看到，CLAHE 中的 clipLimit 参数，最终被转换为了该值乘以 tileSizeTotal （分块像素数）除以 histSize（每个分块的直方图组数），这个转换是干什么呢？是得到每个分组的平均像素数量，如果灰度比较平均的话，每种级别（对应直方图分组数）的灰度所对应的像素数应该相等，当用该平均值乘以 clipLimit，得到的是超过平均值 clipLimit 倍的像素数，这个值就是裁剪的限制值，对于超过这个值的分组就得裁剪。

3.2、直方图裁剪处理过程

裁剪处理在类 CLAHE_Impl 的 apply 方法里调用 CLAHE_CalcLut_Body 类的函数对象来实现的，CLAHE_Impl 是 createCLAHE 生成 CLAHE 实例时真正使用的类，而 CLAHE_CalcLut_Body 类是生成真正的直方图灰度映射和进行裁剪的类。涉及裁剪的代码在 CLAHE_CalcLut_Body 的 operator 函数中。相关的关键源代码如下：

    template <class T, int histSize, int shift>
    void CLAHE_CalcLut_Body<T,histSize,shift>::operator ()(const cv::Range& range) const
    {
			...
            // clip histogram

            if (clipLimit_ > 0)
            {
                // how many pixels were clipped
                int clipped = 0;
                for (int i = 0; i < histSize; ++i)
                {
                    if (tileHist[i] > clipLimit_)
                    {
                        clipped += tileHist[i] - clipLimit_;
                        tileHist[i] = clipLimit_;
                    }
                }

                // redistribute clipped pixels
                int redistBatch = clipped / histSize;
                int residual = clipped - redistBatch * histSize;

                for (int i = 0; i < histSize; ++i)
                    tileHist[i] += redistBatch;

                if (residual != 0)
                {
                    int residualStep = MAX(histSize / residual, 1);
                    for (int i = 0; i < histSize && residual > 0; i += residualStep, residual--)
                        tileHist[i]++;
                }
            }

     		...
        }
    }

从上述代码可以看出，具体裁剪处理时，对于直方图分组像素数超过 clipLimit_的，则将该分组中超过 clipLimit_的像素数累加到 clipped 局部变量中，然后将该直方图分组像素数强制设置为 clipLimit_。

上述过程对当前块的所有分组都处理完成后，将超出后累加的 clipped 变量值按分组数平均分配到各分组中，如果存在不够平均分配的部分，则等间距按顺序插入到分组中，直到所有超出部分都分配到了对应分组。如直方图分组是 256 个，累加的 clipped 值是 1027 个，则先每个分组值加 4，然后将剩余 3 个分别累加到第 0、85、 170 三个位置的直方图分组中。

四、插值处理

插值处理是 CLAHE 算法中理解最困难的，占了本人研究该算法最多的时间，整体算法近 4 个月研究中，插值算法的理解用了 110 多天，也是本人直方图处理一直未能学习完成的根本原因。为了介绍清楚插值处理的算法，下面分成几部分来介绍。

为了说清楚问题，会用到一幅进行直方图均衡处理的经典图像，这幅图像的源图（在老猿的机器上文件名为 f:\pic\valley.png）如下：

4.1、为什么需要插值？

下面这幅图是 clipLimit 设置为 4、tileSize 设置为 4×4（横向和纵向都划分为 4 块，共计 16 块）进行自适应直方图均衡 CLAHE 处理时没有进行插值后的效果图(这是老猿修改算法代码模拟的，因此已经不能称为 CLAHE）：

可以清晰地看到在各个分块之间有明显的灰度突变，整体图像成了棋盘分块一样，这就是图像处理中的“棋盘效应”。为了解决棋盘效应，在 CLAHE 算法中，必须对图像进行插值处理，插值的目的是为了消除各分块之间的突变。下图就是真正的 CLAHE 处理（clipLimit 设置为 4、tileSize 设置为 4×4）后的效果图；

4.2、关于 CLAHE 使用的插值算法

CLAHE 插值处理算法使用的是双线性插值，关于图像处理中的双线性插值网上有一堆的资料介绍，大家可自己查阅，老猿推荐大家阅读《转载：一文讲解图像插值算法原理》。在这些介绍的图像处理中的双线性插值，都是基于四个点的灰度值来计算这个四个点中间的某点的灰度值，但 CLAHE 算法中的双线性插值并不是这样的，这也是老猿理解该算法这么困难的原因。

4.3、CLAHE 双线性插值算法详解

4.3.1、CLAHE 插值算法概述以及插值关联分块

在 CLAHE 算法中进行插值处理是在累积直方图计算和对比度裁剪之后完成的，此时每个分块的灰度映射表已经求出，需要根据每个分块的灰度映射表生成输出图像每个像素的灰度值，生成时使用了双线性插值，也就是要找到与该像素对应的 4 个灰度值以及对应的比率来进行双线性插值。

在具体查找 4 个参考的灰度值以及对应的比例时，CLAHE 是根据当前像素的横坐标 m 以及纵坐标 n 来确认与该像素相关的 4 个分块（这 4 个分块老猿称之为插值关联分块），然后根据（m,n）在源图像的灰度值分别求得在这 4 个分块中直方图均衡后的灰度映射置，再根据该像素横坐标 m 以及点的位置距分块的横向距离在横向进行两次 x 方向的线性插值，再将两次线性插值的结果根据纵坐标 n 及点的位置距分块的纵向距离在纵向进行一次线性插值。

4.3.2、像素 4 个插值关联分块的确定过程

从上面介绍可以知道，需要找到每个输入图像像素点的 4 个分块进行插值，确认这 4 个分块的过程如下：

1. 将图像任意像素点 A（m,n）横坐标减去半个分块宽度、纵坐标减去半个分块高度，得到映射像素点 B（i，j）(注意 B 可能位于图像范围之外，此时会假设图像外存在虚拟扩展的同样大小的分块)。如图为一 4×4 的图像分块中某个像素点 A 和映射点 B(注意一个分块占用 4 个格子，每个格子为 0.5*0.5 分块)：

2. 假设 B 所在分块为 P，其对应分块坐标为(x,y)（对应上图中蓝色标记四个格子）

   • P 右边分块 PR 对应坐标为（x+1,y）、其下边分块 PD 对应坐标(x,y+1)、其斜对角分块 PRD 对应坐标为（x+1,y+1），包括 P 在内这四个分块就是 A 的 4 个插值关联分块；
   • 4 个插值关联分块就是由 2 个横坐标值 x 和 x+1 以及 2 个纵坐标值 y 和 y+1 组合指定的分块；
   • 四个插值关联分块中一定有一个为 A 所在分块、当 B 在图像范围内时一定也包含 B 所在分块，二者可能是同一分块，有可能是不同分块；
   • 上图中这四个图像分块 P(x,y)、PR（x+1,y）、PD(x,y+1)、PRD（x+1,y+1）的分块坐标分别为(1,1)、(1,2)、(2,1)、(2,2)；

3. 当 B 位于图像范围之外时，显然不可能在图像的实际分块范围之内，此时假设在图像范围外也存在虚拟扩展的同样大小的分块

   • 该虚拟扩展分块的分块坐标的 x 值或 y 值至少有一个小于 0；
   • 该虚拟扩展分块并没有真正的图像数据，没有对应的直方图均衡对应的灰度映射表，因此该虚拟扩展分块不能作为插值关联分块；
   • 将该虚拟扩展分块小于 0 的坐标值强制设置为 0，其他非 0 坐标值保持不变，得到的分块坐标对应分块就是插值关联分块。例如当 x=-1 时，关联的分块横坐标 x 和 x+1 值都是 0，假设 y 大于等于 0，则四个插值关联分块及其坐标为 P(0,y)、PR（0,y）、PD(0,y+1)、PRD（0,y+1）；
   • 当 B 位于图像范围之外时，4 个插值关联分块就会包含重复的分块，实际上可能 4 个分块对应的真正分块是 2 个分块（x、y 中只有一个等于-1 时）或者 1 个分块（当 x=y=-1 时）；
   • 这种情况意味着 A 所在范围为距图像边界在半个分块范围内，其目标图像的灰度值不用受其余分块的影响。

下图是一个分块为 4×4 时各个分块中的像素点对应的 4 个插值关联分块的横坐标和纵坐标取值：

注意：

1. 4×4 个分块，在上表中用了 8×8 个格子表示，每个格子表示 0.5 个分块，这样如分块的横坐标为 1，纵坐标为 2 的分块，实际上对应的是第 10 个分块（2*4+2）的 4 个格子，如表格中蓝色底色的 4 格；
2. 上表格子中的 x:(m,n)和 y:(i,j)表示这个格子中的像素的插值关联分块的横坐标为 m 或 n，纵坐标为 i 或 j，对应的插值关联分块为（m，i）、（m，j）、（n，i）、（n，j）4 个分块，当然在靠近边界半格范围内的像素这 4 个分块中有重复分块；
3. 从表格中的数据可以看到，每个分块的像素，当其横坐标在该分块中位于前半部分时，其横向关联的分块是其左边的分块以及其自身所在的分块，当其横坐标在该分块中位于后半部分时，其横向关联的插值关联分块是其右边的分块以及其自身所在的分块；纵坐标类似，只是将关联分块由横向改为了纵向。特别地，对于边界的分块，其上下左右 4 个方向至少有一个方向没有分块，对于这种情况强制设置为该像素所在分块自身；
4. 上表格中的 6 个箭头表示 6 种典型的可能的映射情况，箭头起点代表像素 A（m，n)，终点代表映射点 B(i,j)，按序号 1-6 分别表示：A 和 B 在同一分块、A 和 B 所在分块横坐标相同纵坐标不同、A 和 B 所在分块纵坐标相同横坐标不同、A 和 B 所在分块横坐标纵坐标都不同、B 纵坐标超出图像范围、B 横坐标超出图像范围，当然还有其他情况，但这 6 种情况已经代表了典型的映射情况

4.3.3、线性插值的比例确认

双线性插值实际上包含三次线性插值过程，前两次是 x 方向，最后一次是 y 方向，确认了 4 个插值关联分块之后，分别取原像素的灰度值在这 4 个分块的直方图均衡后灰度映射值（假设为 r00、r01、r10、r11，其中 r00 表示映射点 B 所在分块的直方图均衡后的灰度映射值，r00、r01 对应分块相邻在同一行，r10、r11 对应分块相邻在另一行，而 r00、r10 对应分块相邻在同一列，r01、r11 对应分块相邻在另一列），进行插值时，对于坐标为（m，n）的像素在 x 方向和 y 方向的插值比例是这样确认的：

1. 取映射点 B(i,j)到其所在分块的左边界和右边界与这个分块宽度的比例，作为随后 x 方向的插值比例，同理计算出 y 方向的插值比例；
2. 当(i,j)落在图像外时，此时假设其边界外也存在分块，不过对应方向的分块坐标为-1，其计算插值比例的方法不变。

下面是计算 x 方向插值比例的源代码：

	    float inv_tw = 1.0f / tileSize_.width;
            for (int x = 0; x < src.cols; ++x)
            {
                float txf = x * inv_tw - 0.5f;
                int tx1 = cvFloor(txf);
                 int tx2 = tx1 + 1;
                xa_p[x] = txf - tx1;
                xa1_p[x] = 1.0f - xa_p[x];
                tx1 = std::max(tx1, 0);          
                tx2 = std::min(tx2, tilesX_ - 1);  
            }

上面计算时，x 为输入图像像素的横坐标，tileSize_是分块的大小，tileSize_的 width 为宽度、height 为高度，inv_tw 为分块宽度的倒数，xa_p[x]为映射点到其所在分块左边界的距离占分块宽度的比例，xa1_p[x]为到右边界的距离占分块宽度的比例，tx1 和 tx2 为输入图像像素插值关联分块的 2 个分块横坐标。

下面是计算 y 方向插值比例的源代码：

    float inv_th = 1.0f / tileSize_.height;
    for (int y = range.start; y < range.end; ++y)
    {
        const uchar* srcRow = src_.ptr<uchar>(y);//指向输入图像第y行数据
        uchar* dstRow = dst_.ptr<uchar>(y);//指向输出图像第y行数据
        float tyf = y * inv_th - 0.5f;
        int ty1 = cvFloor(tyf);
        int ty2 = ty1 + 1;
        // 差值作为插值的比例
        float ya = tyf - ty1, ya1 = 1.0f - ya;
        ty1 = std::max(ty1, 0);
        ty2 = std::min(ty2, tilesY_ - 1);
     }

上面计算时，y 为输入图像像素的纵坐标，range.start=0，range.end 等于图像的高度，inv_th 为分块高度的倒数，ty1 和 ty2 为输入图像像素插值关联分块的 2 个分块纵坐标。其中，ya 和 ya1 是映射点到其所在分块上边界和下边界与分块高度的比例。

4.3.4、插值计算

假设输入图像中的像素 A 的映射点到 B 所在分块的分块坐标为(x，y），像素 A 的灰度值在 4 个插值关联分块的灰度映射值为 r00、r01、r10、r11，其位置关系如前所述或编号对应坐标所示。对于 B 所在的分块，点 B 到分块（含虚拟扩展分块）四周的距离与分块宽或高的比例分别是 x1、 x2、y1、y2，显然 x1+x2=1、y1+y2=1。如图：

具体的插值计算过程如下：

1. 进行第一次 x 方向插值，得到 ix1：

ix1 = r00x1+r01x2

2. 进行第二次 x 方向插值，得到 ix2：

ix2 = r10x1+r11x2

3. 进行 y 方向插值，得到输出图像对应的像素的灰度值为 dr：

dr = ix1y1+ix2y2

即：dr = （ r00x1+r01x2）y1+（ r10x1+r11*x2）*y2 。

4.3.4、插值计算过程小结

1. 对输入图像的像素点 A，将其横坐标和纵坐标各减半个分块宽度或半个分块高度，得到映射点 B（映射点 B 可能在图像之外）；

2. 获得 B 距所在分块边界的比例 x1、x2、y1、y2；

3. 根据 B 所在分块的四个边界的 2 个分块横坐标和 2 个分块纵坐标，2 个分块横坐标和 2 个分块纵坐标两两组合得到四个插值关联分块，将 4 个分块中坐标值为-1 的强制设置为 0，此时得到的 4 个坐标对应的 4 个分块为插值关联分块。这 4 个分块在边界处可能重合，此时实际对应 1 个分块（左上角的 1/4 分块）或 2 个分块（其他距离左边界 1/2 宽度内或或上边界的 1/2 高度内的像素）；

4. 根据输入图像像素的灰度值，取该灰度值在四个插值关联分块的各分块直方图均衡灰度映射值，得到 4 个灰度值 r00、r01、r10、r11；

5. 根据公式：dr = （ r00x1+r01x2）y1+（ r10x1+r11*x2）*y2 计算得到输出图像对应像素的灰度值。特别地：

   • 对于 4 个插值关联分块是同一个分块时，r00=r01=r10=r11，x1+x2=1，y1+y2=1，dr 的值就是在该分块内的灰度映射值；

   • 当 4 个插值关联分块是分块纵坐标=0 的分块时（即第一行分块），实际对应两个分块，r00=r10，r01=r11，此时：

     dr =（ r00x1+r01x2）y1+（ r00x1+r01*x2）*y2

     = r00x1(y1+y2)+ r01x2(y1+y2)

     = r00x1+ r01x2

     实际上就是一次 x 方向的插值；

   • 4 个插值关联分块是分块横坐标=0 的分块时（即第一列分块），实际对应两个分块，r00=r01，r10=r11，x1+x2=1，此时：

     dr =（ r00x1+r00x2）y1+（ r10x1+r10*x2）*y2

     = r00y1+r10y2

     实际上就是一次 y 方向的插值。

4.3.5、插值关联分块的获取方法背后的考量

在通过代码解读研究清楚了插值计算过程后，老猿反过头来理解这个计算过程的背后根由，其实这个背后根由很简单：

1. 根据上面给出的棋盘效应的效果图，可以看到将图像分块后各自进行直方图均衡，虽然避免了存在亮区域范围和暗区域范围灰度值变化比较大的情况时整幅图像直方图均衡效果不好的问题，但处理后各分块的灰度不平滑，因此必须进行平滑处理；
2. 当进行插值处理时，要解决的是每个分块之间的灰度值突变问题，而每个分块是进行了直方图均衡处理了的，因此只要确保两个分块的边界的灰度是均衡的，就能保证图像的整体灰度均衡；
3. 假设一个分块其上下左右有 4 个分块，那么将分块等分成 4 个子块，保证左边的两个子块和左边的分块平滑，右边的 2 个子块和右边分块平滑，上边的 2 个子块和上边分块平滑，下边的 2 个子块和下边分块平滑，这样每个子块就会和 2 个相邻的分块平滑，如果同时达到 4 个子块之间平滑，则就能保证该分块和周边 4 个分块平滑；
4. 当对分块 PA 中的像素 A 的横坐标左移半个分块宽度、纵坐标上移半个分块高度得到映射点 B，不考虑 B 在图像外的情况，A 和 B 所在的分块有如下 4 种箭头所示的关系：

• 可以看到，这种移位后的映射点所在分块、左边分块、下边分块和右下对角分块这 4 个插值关联分块是与 A 所在位置最近的几个分块，并且包含了 A 所在分块自身；
• 最妙的是当 B 位于图像之外时，通过强制将小于 0 的分块坐标值设置为 0 后，插值关联分块中就存在重复分块，通过插值关联分块的双线性插值，对这些特殊情况就变成了非插值或者线性插值，有机的将所有情况的插值算法进行了统一，而不用区分所在位置。

五、小结

本文介绍了 OpenCV 对比度受限的自适应直方图均衡 CLAHE 算法，并通过代码研读详细剖了算法中的图像分块整除问题、直方图裁剪过程以及图像插值算法。相关内容在网上公开资料中没有公布，是老猿断断续续花费了 4 个月时间研究剖析，并在春节期间花了大量精力整理的，对计算机视觉算法感兴趣但不熟悉自适应直方图均衡 CLAHE 算法的各位同好应该是非常有帮助的资料。

文章来源：https://xie.infoq.cn/article/f84e1cfee56f9fd6fbe9eb451