自适應直方圖均衡（CLAHE) 代碼及詳細注釋【OpenCV】

效果圖

CLAHE簡介

HE 直方圖增強，大家都不陌生，是一種比較古老的對比度增強算法，它有兩種變體：AHE 和 CLAHE；兩者都是自适應的增強算法，功能差不多，但是前者有一個很大的缺陷，就是有時候會過度放大圖像中相同區域的噪聲，為了解決這一問題，出現了 HE 的另一種改進算法，就是 CLAHE；CLAHE 是另外一種直方圖均衡算法，CLAHE 和 AHE 的差別在于前者對區域對比度實行了限制，并且利用插值來加快計算。它能有效的增強或改善圖像（局部）對比度，進而擷取更多圖像相關邊緣資訊有利于分割。還能夠有效改善 AHE 中放大噪聲的問題。另外，CLAHE 的有一個用途是被用來對圖像去霧。

詳細理論請參考部落格

OpenCV源碼的本地路徑： %OPENCV%\opencv\sources\modules\imgproc\src\clahe.cpp

clahe.cpp

// ----------------------------------------------------------------------
// CLAHE
 
namespace
{
    class CLAHE_CalcLut_Body : public cv::ParallelLoopBody
    {
    public:
        CLAHE_CalcLut_Body(const cv::Mat& src, cv::Mat& lut, cv::Size tileSize, int tilesX, int clipLimit, float lutScale) :
            src_(src), lut_(lut), tileSize_(tileSize), tilesX_(tilesX), clipLimit_(clipLimit), lutScale_(lutScale)
        {
        }
 
        void operator ()(const cv::Range& range) const;
 
    private:
        cv::Mat src_;
        mutable cv::Mat lut_;
 
        cv::Size tileSize_;
        int tilesX_;
        int clipLimit_;
        float lutScale_;
    };
 
    // 計算直方圖查找表
    void CLAHE_CalcLut_Body::operator ()(const cv::Range& range) const
    {
        const int histSize = 256;
 
        uchar* tileLut = lut_.ptr(range.start);
        const size_t lut_step = lut_.step;  // size = tilesX_*tilesY_ * lut_step
 
        // Range(0, tilesX_ * tilesY_)，全圖被分為tilesX_*tiles_Y個塊
        for (int k = range.start; k < range.end; ++k, tileLut += lut_step)
        {
            // (tx, ty)表示目前所在是哪一塊
            // (0, 0) (1, 0)...(tilesX_-1, 0) 
            // (0, 1) (1, 1)...(tilesX_-1, 1) 
            // ...
            // (0, tilesY_-1)... (tilesX_-1, tilesY_-1)
            const int ty = k / tilesX_;
            const int tx = k % tilesX_;
 
            // retrieve tile submatrix
            // 注意：tileSize.width表示分塊的寬度，tileSize.height表示分塊高度
            cv::Rect tileROI;
            tileROI.x = tx * tileSize_.width;   // 換算為全局坐标
            tileROI.y = ty * tileSize_.height;
            tileROI.width = tileSize_.width;
            tileROI.height = tileSize_.height;
 
            const cv::Mat tile = src_(tileROI);
 
            // calc histogram
            int tileHist[histSize] = { 0, };
            // 統計 ROI 的直方圖
            int height = tileROI.height;
            const size_t sstep = tile.step;
            for (const uchar* ptr = tile.ptr<uchar>(0); height--; ptr += sstep)
            {
                int x = 0;
                for (; x <= tileROI.width - 4; x += 4)
                {
                    int t0 = ptr[x], t1 = ptr[x + 1];
                    tileHist[t0]++; tileHist[t1]++;
                    t0 = ptr[x + 2]; t1 = ptr[x + 3];
                    tileHist[t0]++; tileHist[t1]++;
                }
 
                for (; x < tileROI.width; ++x)
                    tileHist[ptr[x]]++;
            }
 
            // 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;
 
                // 處理內插補點的餘數
                for (int i = 0; i < residual; ++i)
                    tileHist[i]++;
            }
 
            // calc Lut
            int sum = 0;
            for (int i = 0; i < histSize; ++i)
            {
                // 累加直方圖
                sum += tileHist[i];
                tileLut[i] = cv::saturate_cast<uchar>(sum * lutScale_); // static_cast<float>(histSize - 1) / tileSizeTotal
            }
        }
    }
 
    class CLAHE_Interpolation_Body : public cv::ParallelLoopBody
    {
    public:
        CLAHE_Interpolation_Body(const cv::Mat& src, cv::Mat& dst, const cv::Mat& lut, cv::Size tileSize, int tilesX, int tilesY) :
            src_(src), dst_(dst), lut_(lut), tileSize_(tileSize), tilesX_(tilesX), tilesY_(tilesY)
        {
        }
 
        void operator ()(const cv::Range& range) const;
 
    private:
        cv::Mat src_;
        mutable cv::Mat dst_;
        cv::Mat lut_;
 
        cv::Size tileSize_;
        int tilesX_;
        int tilesY_;
    };
 
    // 根據相鄰4塊的直方圖插值
    void CLAHE_Interpolation_Body::operator ()(const cv::Range& range) const
    {
        const size_t lut_step = lut_.step;
        // Range(0, src.rows)
        for (int y = range.start; y < range.end; ++y)
        {
            const uchar* srcRow = src_.ptr<uchar>(y);
            uchar* dstRow = dst_.ptr<uchar>(y);
 
            const float tyf = (static_cast<float>(y) / tileSize_.height) - 0.5f;
 
            int ty1 = cvFloor(tyf);
            int ty2 = ty1 + 1;
 
            // 內插補點作為插值的比例
            const float ya = tyf - ty1;
 
            ty1 = std::max(ty1, 0);
            ty2 = std::min(ty2, tilesY_ - 1);
 
            const uchar* lutPlane1 = lut_.ptr(ty1 * tilesX_);   // 目前塊的直方圖
            const uchar* lutPlane2 = lut_.ptr(ty2 * tilesX_);   // 向下一塊的直方圖
 
            for (int x = 0; x < src_.cols; ++x)
            {
                const float txf = (static_cast<float>(x) / tileSize_.width) - 0.5f;
 
                int tx1 = cvFloor(txf);
                int tx2 = tx1 + 1;
 
                // 內插補點作為插值的比例
                const float xa = txf - tx1;
 
                tx1 = std::max(tx1, 0);
                tx2 = std::min(tx2, tilesX_ - 1);
 
                // src_.ptr<uchar>(y)[x]
                const int srcVal = srcRow[x];
 
                // 索引 LUT
                const size_t ind1 = tx1 * lut_step + srcVal;
                const size_t ind2 = tx2 * lut_step + srcVal;  // 向右一塊的直方圖
 
                float res = 0;
                // 根據直方圖的值進行插值
                // lut_.ptr(ty1 * tilesX_)[tx1 * lut_step + srcVa] => lut_[ty1][tx1][srcVal]
                res += lutPlane1[ind1] * ((1.0f - xa) * (1.0f - ya));
                res += lutPlane1[ind2] * ((xa) * (1.0f - ya));
                res += lutPlane2[ind1] * ((1.0f - xa) * (ya));
                res += lutPlane2[ind2] * ((xa) * (ya));
 
                dstRow[x] = cv::saturate_cast<uchar>(res);
            }
        }
    }
 
    class CLAHE_Impl : public cv::CLAHE
    {
    public:
        CLAHE_Impl(double clipLimit = 40.0, int tilesX = 8, int tilesY = 8);
 
        cv::AlgorithmInfo* info() const; // Algorithm類工廠方法封裝相關
 
        void apply(cv::InputArray src, cv::OutputArray dst);
 
        void setClipLimit(double clipLimit);
        double getClipLimit() const;
 
        void setTilesGridSize(cv::Size tileGridSize);
        cv::Size getTilesGridSize() const;
 
        void collectGarbage();
 
    private:
        double clipLimit_;
        int tilesX_;
        int tilesY_;
 
        cv::Mat srcExt_;
        cv::Mat lut_;
    };
 
    CLAHE_Impl::CLAHE_Impl(double clipLimit, int tilesX, int tilesY) :
        clipLimit_(clipLimit), tilesX_(tilesX), tilesY_(tilesY)
    {
    }
        // Algorithm類工廠方法封裝相關 
    //CV_INIT_ALGORITHM(CLAHE_Impl, "CLAHE",
    //  obj.info()->addParam(obj, "clipLimit", obj.clipLimit_);
    //obj.info()->addParam(obj, "tilesX", obj.tilesX_);
    //obj.info()->addParam(obj, "tilesY", obj.tilesY_))
 
        void CLAHE_Impl::apply(cv::InputArray _src, cv::OutputArray _dst)
    {
        cv::Mat src = _src.getMat();
 
        CV_Assert(src.type() == CV_8UC1);
 
        _dst.create(src.size(), src.type());
        cv::Mat dst = _dst.getMat();
 
        const int histSize = 256;
        // 準備 LUT，tilesX_*tilesY_個塊，每個塊都有256個柱子的直方圖
        lut_.create(tilesX_ * tilesY_, histSize, CV_8UC1);
 
        cv::Size tileSize;
        cv::Mat srcForLut;
 
        // 如果分塊剛好（整除）
        if (src.cols % tilesX_ == 0 && src.rows % tilesY_ == 0)
        {
            tileSize = cv::Size(src.cols / tilesX_, src.rows / tilesY_);
            srcForLut = src;
        }
        // 否則對原圖進行擴充
        else
        {
            cv::copyMakeBorder(src, srcExt_, 0, tilesY_ - (src.rows % tilesY_), 0, tilesX_ - (src.cols % tilesX_), cv::BORDER_REFLECT_101);
 
            tileSize = cv::Size(srcExt_.cols / tilesX_, srcExt_.rows / tilesY_);
            srcForLut = srcExt_;
        }
 
        const int tileSizeTotal = tileSize.area();
        const float lutScale = static_cast<float>(histSize - 1) / tileSizeTotal;        // △
 
        // 計算實際的clipLimit
        int clipLimit = 0;
        if (clipLimit_ > 0.0)
        {
            clipLimit = static_cast<int>(clipLimit_ * tileSizeTotal / histSize);
            clipLimit = std::max(clipLimit, 1);
        }
 
        // 分塊并行計算： LUT
        CLAHE_CalcLut_Body calcLutBody(srcForLut, lut_, tileSize, tilesX_, clipLimit, lutScale);
        cv::parallel_for_(cv::Range(0, tilesX_ * tilesY_), calcLutBody);
 
        // 分塊并行計算： 根據直方圖插值
        CLAHE_Interpolation_Body interpolationBody(src, dst, lut_, tileSize, tilesX_, tilesY_);
        cv::parallel_for_(cv::Range(0, src.rows), interpolationBody);
    }
 
    void CLAHE_Impl::setClipLimit(double clipLimit)
    {
        clipLimit_ = clipLimit;
    }
 
    double CLAHE_Impl::getClipLimit() const
    {
        return clipLimit_;
    }
 
    void CLAHE_Impl::setTilesGridSize(cv::Size tileGridSize)
    {
        tilesX_ = tileGridSize.width;
        tilesY_ = tileGridSize.height;
    }
 
    cv::Size CLAHE_Impl::getTilesGridSize() const
    {
        return cv::Size(tilesX_, tilesY_);
    }
 
    void CLAHE_Impl::collectGarbage()
    {
        srcExt_.release();
        lut_.release();
    }
}
 
cv::Ptr<cv::CLAHE> cv::createCLAHE(double clipLimit, cv::Size tileGridSize)
{
    return new CLAHE_Impl(clipLimit, tileGridSize.width, tileGridSize.height);
}

main.cpp

int main(int argc, char** argv)
{
    cv::Mat inp_img = cv::imread("D:/Pictures/beard.jpg");
    if (!inp_img.data) {
        cout << "Something Wrong";
        return -1;
    }
    namedWindow("Input Image", CV_WINDOW_AUTOSIZE);
    cv::imshow("Input Image", inp_img);
 
    cv::Mat clahe_img;
    cv::cvtColor(inp_img, clahe_img, CV_BGR2Lab);
    std::vector<cv::Mat> channels(3);
    cv::split(clahe_img, channels);
 
    cv::Ptr<cv::CLAHE> clahe = cv::createCLAHE();
    // 直方圖的柱子高度大于計算後的ClipLimit的部分被裁剪掉，然後将其平均配置設定給整張直方圖
    // 進而提升整個圖像
    clahe->setClipLimit(4.);    // (int)(4.*(8*8)/256)
    //clahe->setTilesGridSize(Size(8, 8)); // 将圖像分為8*8塊
    cv::Mat dst;
    clahe->apply(channels[0], dst);
    dst.copyTo(channels[0]);
    cv::merge(channels, clahe_img);
 
    cv::Mat image_clahe;
    cv::cvtColor(clahe_img, image_clahe, CV_Lab2BGR);
 
    //cout << cvFloor(-1.5) << endl;
 
    namedWindow("CLAHE Image", CV_WINDOW_AUTOSIZE);
    cv::imshow("CLAHE Image", image_clahe);
    imwrite("out.jpg", image_clahe);
    cv::waitKey(0);
    destroyAllWindows();
    
    return 0;
}

注意：cv::ParallelLoopBody 位于 %OpenCV%\opencv\sources\modules\core\src\parallel.cpp

自适應直方圖均衡（CLAHE) 代碼及詳細注釋【OpenCV】

效果圖

繼續閱讀

查找算法之二分查找查找算法之二分查找

查找算法學習之二分查找（Python版本）——BinarySearch

CQ V1.0分詞bates(基于雙數組tire樹)—應該是目前最快的中文分詞算法

Command Network(POJ 3164)---定根最小樹形圖模闆題題目描述輸入格式輸出格式輸入樣例輸出樣例分析源程式

開源低帶寬語音編解碼器

241 Different Ways to Add Parentheses（C代碼版）

【趨高機器視覺】機器視覺技術原了解析及解決方案

CSMA/CD1． CSMA/CD的概述2． CSMA 的工作原理3． CSMA/CD控制規程及特點4． CSMA/CD協定5． CSMA/CD的優點6．結束語

極大似然法(ML)與最大期望法(EM)

cs231n斯坦福基于卷積神經網絡的CV學習筆記（一）KNN和線性分類器/分類器損失/反向傳播一，KNN圖像分類算法二，線性分類器三，線性分類器損失四，反向傳播五，神經網絡

C++ 第十五周報告1--《冒泡法排序》

筆試面試題目：滑動視窗(二)

資料結構與算法（27）——排序（二）

Dijkstra--簡易版（最短路徑）

GitHub連夜封殺！這份阿裡 10W 字内部 Java 字面試手冊到底有多強？

hdu7108哈希