OpenCV CLAHE直方图均衡解析笔记

Posted 2023-03-01 会飞的鱼chelmx

文章目录

• • 1. CLAHE算法描述
  • 2. CLAHE代码解析

1. CLAHE算法描述

$\text{CLAHE}$是一种限制对比度的自适应直方图均衡算法，关于其详细描述，后续整理给出。

2. CLAHE代码解析

void CLAHE_Impl::apply(cv::InputArray _src, cv::OutputArray _dst)

    CV_INSTRUMENT_REGION();

    CV_Assert( _src.type() == CV_8UC1 || _src.type() == CV_16UC1 );

#ifdef HAVE_OPENCL
    bool useOpenCL = cv::ocl::isOpenCLActivated() && _src.isUMat() && _src.dims()<=2 && _src.type() == CV_8UC1;
#endif

    int histSize = _src.type() == CV_8UC1 ? 256 : 65536;        // 直方图灰度级，灰度图为256

    cv::Size tileSize;                                          // 网格尺寸
    cv::_InputArray _srcForLut;                                 // 预处理图像

    if (_src.size().width % tilesX_ == 0 && _src.size().height % tilesY_ == 0)
                                                                // 判断输入图像是否能按预设网格数目等分
    
        tileSize = cv::Size(_src.size().width / tilesX_, _src.size().height / tilesY_);
                                                                // 网格宽度=图像宽度/网格x方向个数，网格高度=图像高度/网格y方向个数
        _srcForLut = _src;                                      // 将输入图像赋值给预处理图像
    
    else
    
#ifdef HAVE_OPENCL
        if(useOpenCL)
        
            cv::copyMakeBorder(_src, usrcExt_, 0, tilesY_ - (_src.size().height % tilesY_), 0, tilesX_ - (_src.size().width % tilesX_), cv::BORDER_REFLECT_101);
            tileSize = cv::Size(usrcExt_.size().width / tilesX_, usrcExt_.size().height / tilesY_);
            _srcForLut = usrcExt_;
        
        else
#endif
        
            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_);
                                                                // 网格宽度=填充图像宽度/网格x方向个数，网格高度=填充图像高度/网格y方向个数
            _srcForLut = srcExt_;                               // 将填充图像赋值给预处理图像
        
    

    const int tileSizeTotal = tileSize.area();                  // 网格面积
    const float lutScale = static_cast<float>(histSize - 1) / tileSizeTotal;
                                                                // 网格灰度级查找表（概率分布）缩放因子，用于归一化概率分布

    int clipLimit = 0;
    if (clipLimit_ > 0.0)
    
        clipLimit = static_cast<int>(clipLimit_ * tileSizeTotal / histSize);
        clipLimit = std::max(clipLimit, 1);
    

#ifdef HAVE_OPENCL
    if (useOpenCL && clahe::calcLut(_srcForLut, ulut_, tilesX_, tilesY_, tileSize, clipLimit, lutScale) )
        if( clahe::transform(_src, _dst, ulut_, tilesX_, tilesY_, tileSize) )
        
            CV_IMPL_ADD(CV_IMPL_OCL);
            return;
        
#endif

    cv::Mat src = _src.getMat();                                // 输入图像矩阵
    _dst.create( src.size(), src.type() );
    cv::Mat dst = _dst.getMat();                                // 输出图像矩阵
    cv::Mat srcForLut = _srcForLut.getMat();                    // 预处理图像矩阵
    lut_.create(tilesX_ * tilesY_, histSize, _src.type());      // 网格灰度级查找表（概率分布）

    cv::Ptr<cv::ParallelLoopBody> calcLutBody;
    if (_src.type() == CV_8UC1)
        calcLutBody = cv::makePtr<CLAHE_CalcLut_Body<uchar, 256, 0> >(srcForLut, lut_, tileSize, tilesX_, clipLimit, lutScale);
                                                                // 创建灰度级（概率分布）查找表循环主体
    else if (_src.type() == CV_16UC1)
        calcLutBody = cv::makePtr<CLAHE_CalcLut_Body<ushort, 65536, 0> >(srcForLut, lut_, tileSize, tilesX_, clipLimit, lutScale);
    else
        CV_Error( CV_StsBadArg, "Unsupported type" );

    cv::parallel_for_(cv::Range(0, tilesX_ * tilesY_), *calcLutBody);
                                                                // 多线程调用calcLutBody计算各网格灰度级查找表（概率分布）

    cv::Ptr<cv::ParallelLoopBody> interpolationBody;
    if (_src.type() == CV_8UC1)
        interpolationBody = cv::makePtr<CLAHE_Interpolation_Body<uchar, 0> >(src, dst, lut_, tileSize, tilesX_, tilesY_);
                                                                // 创建直方图均衡循环主体
    else if (_src.type() == CV_16UC1)
        interpolationBody = cv::makePtr<CLAHE_Interpolation_Body<ushort, 0> >(src, dst, lut_, tileSize, tilesX_, tilesY_);

    cv::parallel_for_(cv::Range(0, src.rows), *interpolationBody);
                                                                // 多线程调用interpolationBody计算直方图均衡结果

template <class T, int histSize, int shift>                     // <uchar, 256, 0>
void CLAHE_CalcLut_Body<T,histSize,shift>::operator ()(const cv::Range& range) const

    T* tileLut = lut_.ptr<T>(range.start);
    const size_t lut_step = lut_.step / sizeof(T);              // 网格灰度级查找表（概率分布）宽度，等价于histSize（256）

    for (int k = range.start; k < range.end; ++k, tileLut += lut_step)
    
        const int ty = k / tilesX_;                             // 网格y方向编号
        const int tx = k % tilesX_;                             // 网格x方向编号

        // retrieve tile submatrix

        cv::Rect tileROI;                                       // 网格感兴区域
        tileROI.x = tx * tileSize_.width;                       // 网格感兴区域左上角x坐标
        tileROI.y = ty * tileSize_.height;                      // 网格感兴区域左上角y坐标
        tileROI.width = tileSize_.width;                        // 网格感兴区域宽度
        tileROI.height = tileSize_.height;                      // 网格感兴区域高度

        const cv::Mat tile = src_(tileROI);                     // 网格感兴区域矩阵

        // calc histogram

        cv::AutoBuffer<int> _tileHist(histSize);                // 网格直方图
        int* tileHist = _tileHist.data();
        std::fill(tileHist, tileHist + histSize, 0);            // 网格直方图清零

        int height = tileROI.height;                            // 网格感兴区域高度
        const size_t sstep = src_.step / sizeof(T);             // ？？？，等价于图像宽度
        /* 统计感兴区域灰度直方图 */
        for (const T* ptr = tile.ptr<T>(0); height--; ptr += sstep)
        
            int x = 0;
            for (; x <= tileROI.width - 4; x += 4)
            
                int t0 = ptr[x], t1 = ptr[x+1];
                tileHist[t0 >> shift]++; tileHist[t1 >> shift]++;
                t0 = ptr[x+2]; t1 = ptr[x+3];
                tileHist[t0 >> shift]++; tileHist[t1 >> shift]++;
            

            for (; x < tileROI.width; ++x)
                tileHist[ptr[x] >> shift]++;
        

        // clip histogram

        if (clipLimit_ > 0)
        
            // how many pixels were clipped
            /* 将灰度直方图最大频度压缩至clipLimit_ */
            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]++;
            
        

        // calc Lut
        /* 计算灰度级查找表（概率分布） */
        int sum = 0;
        for (int i = 0; i < histSize; ++i)
        
            sum += tileHist[i];
            tileLut[i] = cv::saturate_cast<T>(sum * lutScale_);
        
    

CLAHE_Interpolation_Body(const cv::Mat& src, const cv::Mat& dst, const cv::Mat& lut, const cv::Size& tileSize, const int& tilesX, const int& tilesY) :
    src_(src), dst_(dst), lut_(lut), tileSize_(tileSize), tilesX_(tilesX), tilesY_(tilesY)
                                                                // (src：输入图像矩阵, dst：输出图像矩阵, lut_：网格灰度级查找表（概率分布）, tileSize：网格尺寸, tilesX_：网格x方向个数, tilesY_：网格y方向个数)

    buf.allocate(src.cols << 2);                                // 申请输入图像宽度 * 4大小缓冲区
    ind1_p = buf.data();
    ind2_p = ind1_p + src.cols;
    xa_p = (float *)(ind2_p + src.cols);
    xa1_p = xa_p + src.cols;

    int lut_step = static_cast<int>(lut_.step / sizeof(T));     // 网格灰度级查找表（概率分布）宽度，等价于histSize（256）
    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);                                 // 图像第x + 1列所落入的网格x方向编号
        int tx2 = tx1 + 1;                                      // 图像第x + 1列所落入的网格x方向编号 + 1

        xa_p[x] = txf - tx1;                                    // 图像第x + 1列在网格中距离左边界偏移量
        xa1_p[x] = 1.0f - xa_p[x];                              // 图像第x + 1列在网格中距离下边界偏移量

        tx1 = std::max(tx1, 0);                                 // 修正网格x方向编号（不小于0）
        tx2 = std::min(tx2, tilesX_ - 1);                       // 修正网格x方向编号（不大于网格x方向个数 - 1）

        ind1_p[x] = tx1 * lut_step;                             // 第（tx1，0）个网格在网格灰度级查找表（概率分布）的起始索引
        ind2_p[x] = tx2 * lut_step;                             // 第（tx2，0）个网格在网格灰度级查找表（概率分布）的起始索引
    

void CLAHE_Interpolation_Body<T, shift>::operator ()(const cv::Range& range) const
                                                                // <uchar, 0>
                                                                // range：0 - 输入图像矩阵高度
                                                                // (src：输入图像矩阵, dst：输出图像矩阵, lut_：网格灰度级查找表（概率分布）, tileSize：网格尺寸, tilesX_：网格x方向个数, tilesY_：网格y方向个数)

    float inv_th = 1.0f / tileSize_.height;

    for (int y = range.start; y < range.end; ++y)
    
        const T* srcRow = src_.ptr<T>(y);                       // 输入图像第y + 1行首地址
        T* dstRow = dst_.ptr<T>(y);                             // 输出图像第y + 1行首地址

        float tyf = y * inv_th - 0.5f;

        int ty1 = cvFloor(tyf);                                 // 图像第y + 1行所落入的网格y方向编号
        int ty2 = ty1 + 1;                                      // 图像第y + 1行所落入的网格y方向编号 + 1

        float ya = tyf - ty1, ya1 = 1.0f - ya;                  // 图像第y + 1行在网格中距离上边界以及下边界偏移量

        ty1 = std::max(ty1, 0);                                 // 修正网格y方向编号（不小于0）
        ty2 = std::min(ty2, tilesY_ - 1);                       // 修正网格y方向编号（不大于网格y方向个数 - 1）

        const T* lutPlane1 = lut_.ptr<T>(ty1 * tilesX_);        // 第（0，ty1）个网格灰度级查找表首地址
        const T* lutPlane2 = lut_.ptr<T>(ty2 * tilesX_);        // 第（0，ty2）个网格灰度级查找表首地址

        for (int x = 0; x < src_.cols; ++x)
        
            int srcVal = srcRow[x] >> shift;                    // 获取第（x，y）个像素灰度

            int ind1 = ind1_p[x] + srcVal;                      // 获取第（像素所在网格x方向编号，0）个网格中像素灰度索引
            int ind2 = ind2_p[x] + srcVal;                      // 获取第（像素所在网格x方向编号 + 1，0）个网格中像素灰度索引

            float res = (lutPlane1[ind1] * xa1_p[x] + lutPlane1[ind2] * xa_p[x]) * ya1 +
                        (lutPlane2[ind1] * xa1_p[x] + lutPlane2[ind2] * xa_p[x]) * ya;
                                                                // 取像素所在网格、其右侧、其下侧以及其右下侧灰度均衡结果，并按像素点相对于网格边界的距离进行加权

            dstRow[x] = cv::saturate_cast<T>(res) << shift;     // 将像素处理结果赋值到输出图像