图像处理基础(8)：图像的灰度直方图、直方图均衡化、直方图规定化（匹配）
本文主要介绍了灰度直方图相关的处理，包括以下几个方面的内容：
- 利用OpenCV计算图像的灰度直方图，并绘制直方图曲线
- 直方图均衡化的原理及实现
- 直方图规定化（匹配）的原理及实现
图像的灰度直方图
一幅图像由不同灰度值的像素组成，图像中灰度的分布情况是该图像的一个重要特征。图像的灰度直方图就描述了图像中灰度分布情况，能够很直观的展示出图像中各个灰度级所占的多少。
图像的灰度直方图是灰度级的函数，描述的是图像中具有该灰度级的像素的个数：其中，横坐标是灰度级，纵坐标是该灰度级出现的频率。
不过通常会将纵坐标归一化到\([0,1]\)区间内，也就是将灰度级出现的频率（像素个数）除以图像中像素的总数。灰度直方图的计算公式如下：
其中，\(r_k\)是像素的灰度级，\(n_k\)是具有灰度\(r_k\)的像素的个数，\(MN\)是图像中总的像素个数。
OpenCV灰度直方图的计算
直方图的计算是很简单的，无非是遍历图像的像素，统计每个灰度级的个数。在OpenCV中封装了直方图的计算函数calcHist，为了更为通用该函数的参数有些复杂，其声明如下：
void calcHist( const Mat* images, int nimages,
const int* channels, InputArray mask,
OutputArray hist, int dims, const int* histSize,
const float** ranges, bool uniform = true, bool accumulate = false );
该函数能够同时计算多个图像，多个通道，不同灰度范围的灰度直方图.
其参数如下：
- images，输入图像的数组，这些图像要有相同大大小，相同的深度（
CV_8U CV_16U CV_32F）. - nimages ，输入图像的个数
- channels，要计算直方图的通道个数。
- mask，可选的掩码，不使用时可设为空。要和输入图像具有相同的大小，在进行直方图计算的时候，只会统计该掩码不为0的对应像素
- hist，输出的直方图
- dims，直方图的维度
- histSize，直方图每个维度的大小
- ranges，直方图每个维度要统计的灰度级的范围
- uniform，是否进行归一化，默认为true
- accumulate，累积标志，默认值为false。
为了计算的灵活性和通用性，OpenCV的灰度直方图提供了较多的参数，但对于只是简单的计算一幅灰度图的直方图的话，又显得较为累赘。这里对calcHist进行一次封装，能够方便的得到一幅灰度图直方图。
class Histogram1D
{
private:
int histSize[1]; // 项的数量
float hranges[2]; // 统计像素的最大值和最小值
const float* ranges[1];
int channels[1]; // 仅计算一个通道
public:
Histogram1D()
{
// 准备1D直方图的参数
histSize[0] = 256;
hranges[0] = 0.0f;
hranges[1] = 255.0f;
ranges[0] = hranges;
channels[0] = 0;
}
MatND getHistogram(const Mat &image)
{
MatND hist;
// 计算直方图
calcHist(&image ,// 要计算图像的
1, // 只计算一幅图像的直方图
channels, // 通道数量
Mat(), // 不使用掩码
hist, // 存放直方图
1, // 1D直方图
histSize, // 统计的灰度的个数
ranges); // 灰度值的范围
return hist;
}
Mat getHistogramImage(const Mat &image)
{
MatND hist = getHistogram(image);
// 最大值，最小值
double maxVal = 0.0f;
double minVal = 0.0f;
minMaxLoc(hist, &minVal, &maxVal);
//显示直方图的图像
Mat histImg(histSize[0], histSize[0], CV_8U, Scalar(255));
// 设置最高点为nbins的90%
int hpt = static_cast<int>(0.9 * histSize[0]);
//每个条目绘制一条垂直线
for (int h = 0; h < histSize[0]; h++)
{
float binVal = hist.at<float>(h);
int intensity = static_cast<int>(binVal * hpt / maxVal);
// 两点之间绘制一条直线
line(histImg, Point(h, histSize[0]), Point(h, histSize[0] - intensity), Scalar::all(0));
}
return histImg;
}
};
Histogram1D提供了两个方法：getHistogram返回统计直方图的数组，默认计算的灰度范围是[0,255]；getHistogramImage将图像的直方图以线条的形式画出来，并返回包含直方图的图像。测试代码如下：
Histogram1D hist;
Mat histImg;
histImg = hist.getHistogramImage(image);
imshow("Image", image);
imshow("Histogram", histImg);
其结果如下：
直方图均衡化 Histogram Equalization
假如图像的灰度分布不均匀，其灰度分布集中在较窄的范围内，使图像的细节不够清晰，对比度较低。通常采用直方图均衡化及直方图规定化两种变换，使图像的灰度范围拉开或使灰度均匀分布，从而增大反差，使图像细节清晰，以达到增强的目的。
直方图均衡化，对图像进行非线性拉伸，重新分配图像的灰度值，使一定范围内图像的灰度值大致相等。这样，原来直方图中间的峰值部分对比度得到增强，而两侧的谷底部分对比度降低，输出图像的直方图是一个较为平坦的直方图。
均衡化算法
直方图的均衡化实际也是一种灰度的变换过程，将当前的灰度分布通过一个变换函数，变换为范围更宽、灰度分布更均匀的图像。也就是将原图像的直方图修改为在整个灰度区间内大致均匀分布，因此扩大了图像的动态范围，增强图像的对比度。通常均衡化选择的变换函数是灰度的累积概率，直方图均衡化算法的步骤：
- 计算原图像的灰度直方图 \(P(S_k) = \frac{n_k}{n}\)，其中\(n\)为像素总数，\(n_k\)为灰度级\(S_k\)的像素个数
- 计算原始图像的累积直方图 \(CDF(S_k) = \sum\limits^k_{i=0}\frac{n_i}{n}=\sum\limits^k_{i=0}P_s(S_i)\)
- \(D_j = L\cdot CDF(S_i)\)，其中 \(D_j\)是目的图像的像素，\(CDF(S_i)\)是源图像灰度为i的累积分布，L是图像中最大灰度级（灰度图为255）
其代码实现如下：
- 在上面中封装了求灰度直方图的类，这里直接应用该方法得到图像的灰度直方图；
- 将灰度直方图进行归一化，计算灰度的累积概率；
- 创建灰度变化的查找表
- 应用查找表，将原图像变换为灰度均衡的图像
具体代码如下：
void equalization_self(const Mat &src, Mat &dst)
{
Histogram1D hist1D;
MatND hist = hist1D.getHistogram(src);
hist /= (src.rows * src.cols); // 对得到的灰度直方图进行归一化
float cdf[256] = { 0 }; // 灰度的累积概率
Mat lut(1, 256, CV_8U); // 灰度变换的查找表
for (int i = 0; i < 256; i++)
{
// 计算灰度级的累积概率
if (i == 0)
cdf[i] = hist.at<float>(i);
else
cdf[i] = cdf[i - 1] + hist.at<float>(i);
lut.at<uchar>(i) = static_cast<uchar>(255 * cdf[i]); // 创建灰度的查找表
}
LUT(src, lut, dst); // 应用查找表，进行灰度变化，得到均衡化后的图像
}
上面代码只是加深下对均衡化算法流程的理解，实际在OpenCV中也提供了灰度均衡化的函数equalizeHist,该函数的使用很简单，只有两个参数：输入图像，输出图像。下图为，上述代码计算得到的均衡化结果和调用equalizeHist的结果对比
最左边为原图像，中间为OpenCV封装函数的结果，右边为上面代码得到的结果。
直方图规定化
从上面可以看出，直方图的均衡化自动的确定了变换函数，可以很方便的得到变换后的图像，但是在有些应用中这种自动的增强并不是最好的方法。有时候，需要图像具有某一特定的直方图形状（也就是灰度分布），而不是均匀分布的直方图，这时候可以使用直方图规定化。
直方图规定化，也叫做直方图匹配，用于将图像变换为某一特定的灰度分布，也就是其目的的灰度直方图是已知的。这其实和均衡化很类似，均衡化后的灰度直方图也是已知的，是一个均匀分布的直方图；而规定化后的直方图可以随意的指定，也就是在执行规定化操作时，首先要知道变换后的灰度直方图，这样才能确定变换函数。规定化操作能够有目的的增强某个灰度区间，相比于，均衡化操作，规定化多了一个输入，但是其变换后的结果也更灵活。
在理解了上述的均衡化过程后，直方图的规定化也较为简单。可以利用均衡化后的直方图作为一个中间过程，然后求取规定化的变换函数。具体步骤如下：
- 将原始图像的灰度直方图进行均衡化，得到一个变换函数\(s = T(r)\)，其中s是均衡化后的像素，r是原始像素
- 对规定的直方图进行均衡化，得到一个变换函数\(v = G(z)\)，其中v是均衡化后的像素，z是规定化的像素
- 上面都是对同一图像的均衡化，其结果应该是相等的,\(s = v,且 z = G^{-1}(v) = G^{-1}(T(r))\)
通过，均衡化作为中间结果，将得到原始像素\(r\)和\(z\)规定化后像素之间的映射关系。
详解规定化过程
对图像进行直方图规定化操作，原始图像的直方图和以及规定化后的直方图是已知的。假设\(P_r(r)\)表示原始图像的灰度概率密度，\(P_z(z)\)表示规定化图像的灰度概率密度（r和z分别是原始图像的灰度级，规定化后图像的灰度级）。
- 对原始图像进行均衡化操作，则有\(s_k = T(r_k) = L \cdot \sum\limits_{i=0}^{i=k}P_r(r_k)\)
- 对规定化的直方图进行均衡化操作，则\(v_k = G(z_m) = L \cdot \sum\limits_{j=0}^{j=m}P_z(z_m)\)
- 由于是对同一图像的均衡化操作，所以有\(s_k = v_m\)。
- 规定化操作的目的就是找到原始图像的像素\(s_k\)到规定化后图像像素的\(z_k\)之间的一个映射。有了上一步的等式后，可以得到\(s_k = G(z_k)\)，因此要想找到\(s_k\)想对应的\(z_k\)只需要在\(z\)进行迭代，找到使式子\(G(z_m)-s_k\)的绝对值最小即可。
- 上述描述只是理论的推导过程，在实际的计算过程中，不需要做两次的均衡化操作，具体的推导过程如下：$$
\begin{array}{c}
s_k = v_k \ L \cdot \sum\limits_{i=0}^{i=k}P_r(r_k) = L \cdot \sum\limits_{j=0}^{j=m}P_z(z_m) \ \sum\limits_{i=0}^{i=k}P_r(r_k) = \sum\limits_{j=0}^{j=m}P_z(z_m)
\end
