直方图均衡化

本文主要介绍以下几点内容：

• 直方图是什么
• 直方图均衡化是什么以及有何用途
• 直方图均衡化实现

直方图

直方图是统计一幅图像的每个灰度值出现的频率的图表，直方图反映了图像的像素点灰度值分布的情况，是一种全局的图像特征描述。从数学的角度看，直方图本质是一个关于灰度值的概率分布（pdf)。

直方图均衡化

如果说直方图表示图像中的灰度分布情况，那么直方图均衡化(Histogram Equalization)的作用便是调整灰度值的分布情况，这个“调整”就体现在均衡一词上。直方图均衡化通常在图像的灰度值值分布过于集中，即直方图比较陡的情况下。均衡化就是将灰度的分布从集中调整为均匀分布，经此处理好，在图像的表现为：图像高亮部分和低暗部分减少，图像的对比度提高，之前过曝或者曝光不足之处细节都得到提升。所以直方图均衡化对于背景和前景太亮或者太暗的图像是非常有效的。

那直方图均衡化调整的是如何调整灰度分布的呢，其数学本质又是什么。原理其实非常简单，直方图是灰度份概率分布（pdf),而直方图均衡化实际就是将原图像的灰度概率分布图pdf进行变换调整得到一个新的pdf。实现原理为对此概率分布求积分，得到灰度的累计分布，即cdf。

直方图均衡化步骤

算法步骤：

（1） 对于有G个灰度级（一般是256）大小为MxN的图像，创建一个长为G数组H并初始化为0.

（2） 形成图像直方图：扫描每个像素，增加相应的H成员，当像素p具有亮度$g_p$d时，做
$$H[g_p]=H[g_p]+1$$

说明：以上便是求取直方图的运算，注意直方图无需归一化，此处直方图表示的是每个灰度级出现的次数而非概率。
（3）形成累积的直方图$H_c$:
$$H_c[0] = H_c[0]$$
$$H_c[p] = H_c[p-1] + H[p] \quad p = 1,2,...,G-1$$

说明：以上便是将直方图做离散的积分运算得到的累积直方图，即cdf（同样没归一化），此步为均衡化的关键步骤，由公式可知，将$H_c[p]$灰度级出现次数做调整得到一个新的$H_c[p]$。实际上此步骤已经得到了新的灰度值直方图，接下来的步骤是将新直方图应用到图像中。
（4）置$T[p] = round(\frac{G-1}{NM}H_c[p])$。这一步骤构造了一个是MN倍数的与单调增加的$H_c$中的值对应的查找表，有助于提高实现的效率。

说明：此步的输入为一个灰度值$p$, 结果输出为一个新的灰度值。分两步来理解，首先$\frac{H_c[p]}{NM}$为归一化运算，为$p$的出现概率。然后$(G-1)\times\frac{H_c[p]}{NM}$,灰度级乘以概率得到新的灰度值。至此就完成了一个像素点灰度值的调整。

（5）重新扫描图像，写一个具有灰度级$g_p$的输出图像，设置
$$g_p = T[g_p]$$

说明：此步骤为将图像将每个像素点的像素值按照新的直方图做调整。

代码实现

// Copyright https://mangoroom.cn
// License(MIT)
// Author:mango
// Histogram Equalization
// this is HistogramEqualization.cpp

#include <iostream>

#include<opencv2/opencv.hpp>
#include<array>

namespace imageprocess
{
    // gray histogram 
    void GrayHistogram(const cv::Mat& gray_image, std::array<int, 256>& histogram);

    // Histogram equalization
    void HistogramEqualization(const std::array<int, 256>& histogram, std::array<int, 256>& out, int pixels_cout);

    // histogram array to Mat
    void Histogram2Mat(const std::array<int, 256>& histogram, cv::Mat& histogram_mat);

}//namespace imageproccess

int main()
{
    cv::Mat src_image = cv::imread("Fig0222(a)(face).tif", cv::IMREAD_GRAYSCALE);

    if (src_image.empty())
    {
        return -1;
    }

    std::array<int, 256> histogram = { 0 };
    std::array<int, 256> new_histogram = { 0 };
    std::array<int, 256> dst_histogram = { 0 };

    imageprocess::GrayHistogram(src_image, histogram);

    cv::Mat histogram_mat;
    cv::Mat cdf;
    cv::Mat new_histogram_mat;
    cv::Mat dst_image = src_image.clone();
    imageprocess::HistogramEqualization(histogram, new_histogram, src_image.rows * src_image.cols);
    imageprocess::Histogram2Mat(histogram, histogram_mat);
    imageprocess::Histogram2Mat(new_histogram, cdf);

    // adjust the origin image pixels
    for (size_t i = 0; i < src_image.rows; i++)
    {
        for (size_t j = 0; j < src_image.cols; j++)
        {
            dst_image.at<uchar>(i, j) = new_histogram.at(src_image.at<uchar>(i, j));
        }
    }

    imageprocess::GrayHistogram(dst_image, dst_histogram);
    imageprocess::Histogram2Mat(dst_histogram, new_histogram_mat);

    cv::imshow("src-image", src_image);
    cv::imshow("dst-image", dst_image);
    cv::imshow("histogram", histogram_mat);
    cv::imshow("new-histogram", new_histogram_mat);
    cv::imshow("cdf", cdf);

    cv::imwrite("src-image.jpg", src_image);
    cv::imwrite("dst-image.jpg", dst_image);
    cv::imwrite("histogram.jpg", histogram_mat);
    cv::imwrite("new-histogram.jpg", new_histogram_mat);
    cv::imwrite("cdf.jpg", cdf);

    cv::waitKey(0);

    return 0;
}

void imageprocess::HistogramEqualization(const std::array<int, 256>& histogram, std::array<int, 256>& out, int pixels_cout)
{
    // check the input parameter
    assert(!histogram.empty() && !out.empty());

    // calculate the new histogram (cdf)
    out.at(0) = histogram.at(0);
    for (size_t i = 1; i < 256; i++)
    {
        out.at(i) = out.at(i - 1) + histogram.at(i);
    }

    // create the look up table
    for (size_t i = 0; i < 256; i++)
    {
        out.at(i) = static_cast<int>(255.0 * out.at(i) / pixels_cout);
    }
}


void imageprocess::GrayHistogram(const cv::Mat& gray_image, std::array<int, 256>& histogram)
{
    // check the input parameter : 检查输入参数
    assert(gray_image.channels() == 1);
    assert(histogram.size() == 256);

    // step1: All elements of the histogram array are assigned a value of 0 : 将数组histogram所有的元素赋值为0
    histogram = { 0 };

    // step2: Do hf[f(x,y)]+1 for all pixels of the image: 对图像所有元素，做hf[f(x,y)]+1
    for (size_t i = 0; i < gray_image.rows; i++)
    {
        for (size_t j = 0; j < gray_image.cols; j++)
        {
            int z = gray_image.at<uchar>(i, j);
            histogram.at(z) += 1;
        }
    }
}

void imageprocess::Histogram2Mat(const std::array<int, 256>& histogram, cv::Mat& histogram_mat)
{
    // Check the input parameter :检查输入参数
    assert(histogram.size() == 256);

    // step1: calculate the row of mat : 计算mat的row值
    int row = 0;
    for (size_t i = 0; i < histogram.size(); i++)
    {
        row = row > histogram.at(i) ? row : histogram.at(i);
    }

    // step2: initialize mat : 初始化mat
    histogram_mat = cv::Mat::zeros(row, 256, CV_8UC1);

    // step3: assign value for mat : 为mat赋值
    for (size_t i = 0; i < 256; i++)
    {
        int gray_level = histogram.at(i);

        if (gray_level > 0)
        {
            histogram_mat.col(i).rowRange(cv::Range(row - gray_level, row)) = 255;
        }
    }

    // step4: resize the histogram mat : 缩放直方图
    cv::resize(histogram_mat, histogram_mat, cv::Size(256, 256));
}

references

【1】图像处理、分析与机器视觉4th-page81
【2】数字图像处理3rd-page74
【3】直方图均衡化-维基百科

本文由芒果浩明发布，转载请注明来源。
本文链接：https://mangoroom.cn/opencv/histogram-equalization.html