本节为图像空域增强的中篇,介绍一类基础图像空域增强技术:直方图操作。
通过本节的学习,读者将初步了解图像空域增强中常用的直方图操作,比如直方图、归一化直方图和累计直方图的计算方法,以及直方图均衡、直方图匹配、CLAHE等三种典型直方图操作。这一部分主要用到了Skimage中的exposure模块中的部分函数。
目录
5.3 直方图操作
5.3.1 直方图、归一化直方图、累计直方图
图像直方图(Histogram)是反映图像像素分布的统计表,记为ℎ(𝑘),横坐标代表像素值的取值区间,纵坐标代表每一像素值在图像中的像素总数(对应普通直方图)或者所占的百分比(对应归一化直方图)。
归一化直方图(Uniform Histogram)是图像灰度级的函数,用来描述每个灰度级在图像矩阵中的发生的概率,记为𝑝(𝑘),它可以ℎ(𝑘)用得到。
累计分布直方图(Cumulitive Distribution Histogram),简称累计直方图,代表图像组成成分在灰度级的累计概率分布情况,每一个概率值代表小于等于此灰度值的概率,记为𝑐(𝑘),它也可以用ℎ(𝑘)得到。
以下是三种直方图的定义:
式中,𝑘是图像灰度级的取值,表示取值为𝑘的像素点数目,N是图像内像素点的总数。
下图是多幅图像及其对应的直方图情况,不难看出,灰度直方图反映了图像中的灰度分布规律,直观地表现了图像中各灰度级的占比,很好地体现出图像的亮度和对比度信息:
- 灰度图分布居中说明亮度正常,偏左说明亮度较暗,偏右表明亮度较高;
- 狭窄陡峭表明对比度降低,宽泛平缓表明对比度较高。
(从左至右分别是较暗、较亮、较黯淡,以及高对比度图像及其直方图示例,改编自参考文献1)
在Skimage中使用histogram函数计算图像的直方图。
histogram函数完成反色处理。
看一下histogram函数的声明:
exposure.histogram(image, nbins, source_range, normalize, channel_axis)
部分参数说明:
- image:输入图像。
- nbins:用于计算直方图的bin的数目。
- source_range:待补充。
- normalize:待补充 。
- channel_axis:待补充。
返回值:
- hist:灰度直方图的数值,数组类型。如果通道数不止一个,则hist是二维数组。
- bin_centers:每个bin的中心值,数组类型。
在Skimage中使用cumulative_distribution函数计算图像的累计直方图。
cumulative_distribution函数完成反色处理。
看一下cumulative_distribution函数的声明:
exposure.cumulative_distribution(image, nbins)
部分参数说明:
- image:输入图像。
- nbins:用于计算累计直方图的bin的数目,默认值是256。
返回值:
- img_cdf:累计直方图,数组类型。
- bin_centers:每个bin的中心值,数组类型。
以下是一段计算并显示图像直方图的代码(需修订)
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from skimage import data, img_as_float
from skimage import exposure
matplotlib.rcParams['font.size'] = 8
def plot_img_and_hist(image, axes, bins=256):
"""Plot an image along with its histogram and cumulative histogram.
"""
image = img_as_float(image)
ax_img, ax_hist = axes
ax_cdf = ax_hist.twinx()
# Display image
ax_img.imshow(image, cmap=plt.cm.gray)
ax_img.set_axis_off()
# Display histogram
ax_hist.hist(image.ravel(), bins=bins, histtype='step', color='black')
ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
ax_hist.set_xlabel('Pixel intensity')
ax_hist.set_xlim(0, 1)
ax_hist.set_yticks([])
# Display cumulative distribution
img_cdf, bins = exposure.cumulative_distribution(image, bins)
ax_cdf.plot(bins, img_cdf, 'r')
ax_cdf.set_yticks([])
return ax_img, ax_hist, ax_cdf
image = img_as_float(data.camera())
hk = exposure.histogram(image)
cdf = exposure.cumulative_distribution(image)
5.3.2 直方图均衡
直方图均衡化(Histogram Equalization,简称直方图均衡,HE)是一种简单有效的图像增强技术。根据直方图的形态可以判断图像的质量,通过调控直方图的形态可以改善图像的质量。
直方图均衡的基本思想是对图像中占比大的灰度级进行展宽,而对占比小的灰度级进行压缩,使图像的直方图分布较为均匀,扩大灰度值差别的动态范围,从而增强图像整体的对比度。
因此,直方图均衡就是对图像进行非线性拉伸,重新分配图像像素值,本质上是根据直方图对图像进行线性或非线性灰度变换。例如,直方图均衡化可以把原始图像的直方图调整到均匀分布,增加像素之间灰度值差别的动态范围,从而增强图像整体的对比度。
直方图均衡的转换函数需要用到累计直方图,详细推导过程可阅读参考文献1,在此列出主要步骤:
- 计算原始灰度图像的直方图ℎ(𝑘);
- 计算原始图像的累计分布函数𝑐(𝑘);
- 通过插值计算得到新的灰度值𝑘′。
在Skimage中使用equalize_hist函数计算图像的累计直方图。
equalize_hist函数完成反色处理。
看一下equalize_hist函数的声明:
exposure.equalize_hist(image, nbins, mask)
部分参数说明:
- image:输入图像。
- nbins:用于计算累计直方图的bin的数目,默认值是256。
- mask:二值模板图像,和原图大小相同,取值为True的像素点才进行直方图均衡处理。
返回值:
- out:均衡化处理结果。
以下是一段计算并显示图像直方图的代码(需修订)
5.3.3 直方图匹配
直方图均衡直接对图像全局进行均衡化,生成具有均匀直方图的图像,但并没有考虑局部图像区域的具体情况,有时需要生成具有特殊形状直方图的图像。
直方图匹配(Histogram mapping)又称为直方图规定化(Histogram Specification),是指将图像的直方图调整为规定的形状。 例如,将一幅图像或某一区域的直方图匹配到另一幅影像上,使两幅影像的色调保持一致。
直方图匹配是建立在直方图均衡处理的基础上。用𝑟和𝑧分别表示输入图像匹配图像,和
分别表示输入图像和匹配图像的直方图。则基于直方图均衡技术的直方图匹配主要步骤是:
- 通过
,计算图像𝑟的直方图均衡变换函数𝑇,满足
;
- 通过
计算图像z的直方图均衡变换函数𝐺,满足
;
- 计算变换函数𝐺的逆变换函数
,即满足
;
- 对输入图像𝑟进行直方图均衡得到均衡图像𝑠,然后再用逆变换函数
将其映射到
,得到直方图匹配图像𝑧。本步骤中的两次变换,也可以合并为一次完成,即总的变换为
。
有关直方图匹配算法的详细介绍,可查阅参考文献1。
在Skimage中使用match_histograms函数计算图像的累计直方图。
match_histograms函数完成反色处理。
看一下match_histograms函数的声明:
skimage.exposure.match_histograms(image, reference, channel_axis)
部分参数说明:
- image:输入图像。
- reserence:参考图像,其通道数需要如输入图像相同。
- channel_axis:指定用于进行匹配化处理的通道序号,如果不指定,则默认为灰度图像。
返回值:
- out:直方图匹配处理结果。
下面给出一段完成直方图匹配化的代码
import matplotlib.pyplot as plt
from skimage import data
from skimage import exposure
from skimage.exposure import match_histograms
reference = data.coffee()
image = data.chelsea()
matched = match_histograms(image, reference, channel_axis=-1)
fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(8, 3),
sharex=True, sharey=True)
for aa in (ax1, ax2, ax3):
aa.set_axis_off()
ax1.imshow(image)
ax1.set_title('Source')
ax2.imshow(reference)
ax2.set_title('Reference')
ax3.imshow(matched)
ax3.set_title('Matched')
plt.tight_layout()
plt.show()
直方图匹配结果如下,其中从左至右分别是:原图、参考图和匹配结果图。
(Histogram matching — skimage 0.22.0 documentation (scikit-image.org))
5.3.4 对比度受限自适应直方图均衡(CLAHE)
直方图均衡和直方图匹配都是基于整幅图像的灰度分布进行全局变换,并非针对图像局部区域的细节进行增强。一方面,传统的直方图均衡技术由于是利用了整幅图像灰度值的统计信息,因此在变换时没有考虑空间信息,在某些情况下效果一般;另一方面,直方图均衡处理对于局部同样适用,局部直方图处理的思想是基于像素邻域的灰度分布进行直方图变换处理。
局部直方图处理过程是:
- 设定某一大小的模板(矩形邻域),在图像中沿逐个像素移动;
- 对每个像素位置,计算模板区域的直方图,对该局部区域进行直方图均衡或直方图匹配变换,变换结果只用于模板区域中心像素点的灰度值修正;
- 模板(邻域)在图像中逐行逐列移动,遍历所有像素点,完成对整幅图像的局部直方图处理。
Skimage提供了一种经典的局部自适应直方图均衡方法的实现,它称为“限制对比度自适应直方图均衡化”方法(Contrast Limited Adaptive Hitogram Equalization:CLAHE),采用了限制直方图分布的方法和加速的插值方法。有关CLAHE的原理,可查阅参考文献2,在此不再赘述。
在Skimage中使用equalize_adapthist函数实现CLAHE算法。
equalize_adapthist函数完成CLAHE处理。
看一下equalize_adapthist函数的声明:
exposure.equalize_adapthist(image, kernel_size, clip_limit, nbins)
部分参数说明:
- image:输入图像。
- kernel_size:该算法所用模板大小。
- clip_limit:截断参数,取值在[0,1]之间,取值越高,处理结果的对比度越大。
返回值:
- out:自适应直方图均衡化处理结果。
有关直方图均衡和自适应直方图均衡处理结果的对比,参考下一小节。
5.3.5 综合实例
以下给出一段用三种典型灰度变换技术处理低质图像的实例,所用增强技术包括了上一节介绍的点运算中的对比度拉伸,以及本节介绍的直方图均衡和自适应直方图均衡化处理。
# Load an example image
img = data.moon()
# Contrast stretching
p2, p98 = np.percentile(img, (2, 98))
img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))
# Equalization
img_eq = exposure.equalize_hist(img)
# Adaptive Equalization
img_adapteq = exposure.equalize_adapthist(img, clip_limit=0.03)
# Display results
fig = plt.figure(figsize=(8, 5))
axes = np.zeros((2, 4), dtype=object)
axes[0, 0] = fig.add_subplot(2, 4, 1)
for i in range(1, 4):
axes[0, i] = fig.add_subplot(2, 4, 1+i, sharex=axes[0,0], sharey=axes[0,0])
for i in range(0, 4):
axes[1, i] = fig.add_subplot(2, 4, 5+i)
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0])
ax_img.set_title('Low contrast image')
y_min, y_max = ax_hist.get_ylim()
ax_hist.set_ylabel('Number of pixels')
ax_hist.set_yticks(np.linspace(0, y_max, 5))
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_rescale, axes[:, 1])
ax_img.set_title('Contrast stretching')
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_eq, axes[:, 2])
ax_img.set_title('Histogram equalization')
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_adapteq, axes[:, 3])
ax_img.set_title('Adaptive equalization')
ax_cdf.set_ylabel('Fraction of total intensity')
ax_cdf.set_yticks(np.linspace(0, 1, 5))
# prevent overlap of y-axis labels
fig.tight_layout()
plt.show()
(处理结果从左至右依次是:低对比度原图、对比度拉伸、直方图均衡、自适应直方图均衡)
Histogram Equalization — skimage 0.22.0 documentation (scikit-image.org)
参考文献:
【1】《数字图像处理》(第4版)冈萨雷斯等,电子工业出版社,2020年.
【2】Zuiderveld, Karel. "Contrast limited adaptive histogram equalization." Graphics gems IV. Academic Press Professional, Inc., 1994.
(本节初稿完成时间:2024-02-12)
(欢迎对DIP+python算法开发感兴趣的初学者,尤其是相关专业本科和低年级研究生关注,本专栏完将持续更新,总篇数不会少于50篇,每篇不会少于5000字,专栏完成之前(估计至少半年)完全免费阅读,敬请关注)