本文主要介绍了灰度直方图相关的处理,包括以下几个方面的内容:
一幅图像由不同灰度值的像素组成,图像中灰度的分布情况是该图像的一个重要特征。图像的灰度直方图就描述了图像中灰度分布情况,能够很直观的展示出图像中各个灰度级所占的多少。
图像的灰度直方图是灰度级的函数,描述的是图像中具有该灰度级的像素的个数:其中,横坐标是灰度级,纵坐标是该灰度级出现的频率。
不过通常会将纵坐标归一化到[0,1]区间内,也就是将灰度级出现的频率(像素个数)除以图像中像素的总数。灰度直方图的计算公式如下:
p(r_k)=\frac{n_k}{MN}
其中,rk是像素的灰度级,nk是具有灰度rk的像素的个数,MN是图像中总的像素个数。
直方图的计算是很简单的,无非是遍历图像的像素,统计每个灰度级的个数。在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 );
该函数能够同时计算多个图像,多个通道,不同灰度范围的灰度直方图. 其参数如下
为了计算的灵活性和通用性,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);
假如图像的灰度分布不均匀,其灰度分布集中在较窄的范围内,使图像的细节不够清晰,对比度较低。通常采用直方图均衡化及直方图规定化两种变换,使图像的灰度范围拉开或使灰度均匀分布,从而增大反差,使图像细节清晰,以达到增强的目的。
直方图均衡化,对图像进行非线性拉伸,重新分配图像的灰度值,使一定范围内图像的灰度值大致相等。这样,原来直方图中间的峰值部分对比度得到增强,而两侧的谷底部分对比度降低,输出图像的直方图是一个较为平坦的直方图。
直方图的均衡化实际也是一种灰度的变换过程,将当前的灰度分布通过一个变换函数,变换为范围更宽、灰度分布更均匀的图像。也就是将原图像的直方图修改为在整个灰度区间内大致均匀分布,因此扩大了图像的动态范围,增强图像的对比度。通常均衡化选择的变换函数是灰度的累积概率,直方图均衡化算法的步骤:
P(S_k) = \frac{n_k}{n}
,其中n为像素总数,Nk为灰度级Sk的像素个数
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)
其中Dj为目的图像的像素,
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))
对图像进行直方图规定化操作,原始图像的直方图和以及规定化后的直方图是已知的。假设
P_r(r)
表示原始图像的灰度概率密度,Pz(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
\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{array}
上面公式表示,假如
规定化后的对应灰度是
的话,需要满足的条件是
的累积概率和
的累积概率是最接近的。 下面是一个具体计算的例子:
首先得到原直方图的各个灰度级的累积概率\(V_s\)以及规定化后直方图的各个灰度级的累积概率\(V_z\),那么确定\(s_k\)到\(z_m\)之间映射关系的条件就是:
的值最小。 以\(k = 2\)为例,其原始直方图的累积概率是:0.65,在规定化后的直方图的累积概率中和0.65最接近(相等)的是灰度值为5的累积概率密度,则可以得到原始图像中的灰度级2,在规定化后的图像中的灰度级是5。
直方图规定化的实现可以分为一下三步:
具体代码实现如下:
void hist_specify(const Mat &src, const Mat &dst,Mat &result)
{
Histogram1D hist1D;
MatND src_hist = hist1D.getHistogram(src);
MatND dst_hist = hist1D.getHistogram(dst);
float src_cdf[256] = { 0 };
float dst_cdf[256] = { 0 };
// 源图像和目标图像的大小不一样,要将得到的直方图进行归一化处理
src_hist /= (src.rows * src.cols);
dst_hist /= (dst.rows * dst.cols);
// 计算原始直方图和规定直方图的累积概率
for (int i = 0; i < 256; i++)
{
if (i == 0)
{
src_cdf[i] = src_hist.at<float>(i);
dst_cdf[i] = dst_hist.at<float>(i);
}
else
{
src_cdf[i] = src_cdf[i - 1] + src_hist.at<float>(i);
dst_cdf[i] = dst_cdf[i - 1] + dst_hist.at<float>(i);
}
}
// 累积概率的差值
float diff_cdf[256][256];
for (int i = 0; i < 256; i++)
for (int j = 0; j < 256; j++)
diff_cdf[i][j] = fabs(src_cdf[i] - dst_cdf[j]);
// 构建灰度级映射表
Mat lut(1, 256, CV_8U);
for (int i = 0; i < 256; i++)
{
// 查找源灰度级为i的映射灰度
// 和i的累积概率差值最小的规定化灰度
float min = diff_cdf[i][0];
int index = 0;
for (int j = 1; j < 256; j++)
{
if (min > diff_cdf[i][j])
{
min = diff_cdf[i][j];
index = j;
}
}
lut.at<uchar>(i) = static_cast<uchar>(index);
}
// 应用查找表,做直方图规定化
LUT(src, lut, result);
}
上面函数的第二个参数的直方图就是规定化的直方图。代码比较简单,这里就不一一解释了。其结果如下:
左边是原图像,右边是规定化的图像,也就是上面函数的第一个和第二个输入参数。原图像规定化的结果如下:
原图像规定化后的直方图和规定化的图像的直方图的形状比较类似, 并且原图像规定化后整幅图像的特征和规定化的图像也比较类似,例如:原图像床上的被子,明显带有规定化图像中水的波纹特征。
直方图规定化过程中,在做灰度映射的时候,有两种常用的方法:
对于GML的映射方法,一直没有很好的理解,但是根据其算法描述实现了该方法,代码这里先不放出,其处理结果如下:
其结果较SML来说更为亮一些,床上的波浪特征也更为明显,但是其直方图形状,和规定化的直方图对比,第一个峰不是很明显。