天津中医药大学-20级医学信息工程 教师:王翌 学生:邓集亲
声明:本文章中所涉及的代码并非我个人独立完成的成果。
在撰写的过程中,我广泛地吸取了前一辈人,尤其是学长学姐们的宝贵学习经验。通过深入研究他们的学习轨迹,以及查阅和分析了众多相关的理论文献与资料,并在王老师的悉心指导下,经过反复的实验、调试与优化,最终得以总结完成本文所展现的代码。
建议各位学弟学妹们,先根据王老师的授课内容,独立思考一下应该如何完成。如果实在是有理解上的困难,不知道从何下手,再来学习和参考本文。
学长我是用Python写的,如果你使用MATLAB,同样可以参考此代码进行翻译。我在代码中加入了许多注释和测试环节,以便于理解。由于学长能力有限,代码中或许存在一些疏漏或错误,请批判性地参考。
实验三:直方图均衡化
作业要求:
- 参考“空间域图像增强”和“彩色图像处理”课的内容,对灰度和彩色图片进行直方图均衡化处理,输出均衡化后的图片。
- 分别在RGB颜色空间和HSI颜色空间下对彩色图片进行直方图均衡化操作,输出均衡化后的图片,观察并分析两个颜色空间下实验结果的差别。
例子:
例一:直方图均衡化:
例二:HSI颜色空间下的直方图均衡化:
实验图片:
灰度图片直方图均衡化:histeq1.jpg、histeq2.jpg、histeq3.jpg、histeq4.jpg。
彩色图片直方图均衡化:histeqColor.jpg。
再自选其它至少20张图片(包括灰度图片和彩色图片)。
图像的灰度直方图描述了图像中灰度分布情况, 能够很直观的展示出图像中各个灰度级所占的多少。图像的灰度直方图是灰度级的函数, 描述的是图像中具有该灰度级的像素的个数: 其中, 横坐标是灰度级, 纵坐标是该灰度级出现的率。
它描述每个灰度级具有的像素个数, 但不包含这些像素在图像中的位置信息。 图像直方图不关心像素所处的空间位置, 因此不受图像旋转和平移变化的影响, 可以作为图像的特征。任何一幅特定的图像都有唯一的直方图与之对应, 但不同的图像可以有相同的直方图。如果一幅图像有两个不相连的区域组成, 并且每个区域的直方图已知, 则整幅图像的直方图是该两个区域的直方图之和。
直方图均衡化是一种增强图像对比度的方法,其主要思想是将一副图像的直方图分布通过累积分布函数变成近似均匀分布,从而增强图像的对比度。为了将原图像的亮度范围进行扩展, 需要一个映射函数, 将原图像的像素值均衡映射到新直方图中,即累积分布函数。
Windows图片文件目录:
ImageSet文件目录下储存其他待处理图片
Python代码:
#导入库
import cv2
import numpy as np
import matplotlib.pyplot as plt
import collections
#直方图均衡化
def histogram_equalization(histogram_e, lut_e, image_e):
sum_temp = 0
cf = []
for i in histogram_e:
sum_temp += i
cf.append(sum_temp)
for i, v in enumerate(lut_e):
lut_e[i] = int(255.0 * (cf[i] / sum_temp) + 0.5)
equalization_result = lut_e[image_e]
return equalization_result
# 计算灰度图的直方图
def draw_histogram(grayscale):
gray_key = []
gray_count = []
gray_result = []
histogram_gray = list(grayscale.ravel()) # 将多维数组转换成一维数组
gray = dict(collections.Counter(histogram_gray)) # 统计图像中每个灰度级出现的次数
gray = sorted(gray.items(), key=lambda item: item[0]) # 根据灰度级大小排序
for element in gray:
key = list(element)[0]
count = list(element)[1]
gray_key.append(key)
gray_count.append(count)
for i in range(0, 256):
if i in gray_key:
num = gray_key.index(i)
gray_result.append(gray_count[num])
else:
gray_result.append(0)
gray_result = np.array(gray_result)
return gray_result
x = []
for i in range(0, 256): # 横坐标
x.append(i)
#读取图像
original = cv2.imread(r'D:\deng\smalljob\ImageSet\histeq1.jpg')
original_gray = cv2.cvtColor(original, cv2.COLOR_BGR2GRAY) #将图像从RGB颜色空间转换到灰度颜色空间
histogram = draw_histogram(original_gray) # 直方图转化
plt.bar(x, histogram) # 绘制原图直方图
plt.savefig('D:/deng/before_histogram.png')# 这里暂时不关闭绘制直方图的窗口,在处理图片后,再在该窗口绘制直方图作对比
before_histogram=cv2.imread(r'D:\deng\before_histogram.png')
lut = np.zeros(256, dtype=original_gray.dtype) # 创建空的查找表,返回image类型的用0填充的数组;
result = histogram_equalization(histogram, lut, original_gray) # 均衡化处理
histogram_equ = draw_histogram(result) # 直方图转化
plt.bar(x, histogram_equ) # 绘制处理后图像直方图
plt.savefig('D:/deng/after_histogram.png')
plt.close()#关闭绘制直方图的窗口
after_histogram=cv2.imread(r'D:\deng\after_histogram.png')
#展示结果
plt.subplot(221),plt.imshow(original,cmap='gray')
plt.title('original')
plt.axis('off')
plt.subplot(222),plt.imshow(before_histogram,cmap='gray')
plt.title('before_histogram')
plt.axis('off')
plt.subplot(223),plt.imshow(result,cmap='gray')
plt.title('result')
plt.axis('off')
plt.subplot(224),plt.imshow(after_histogram,cmap='gray')
plt.title('after_histogram')
plt.axis('off')
plt.show()
运行结果:
以上是处理灰度图,以下是处理彩色图像(RGB颜色模式下)
#导入库
import cv2
import numpy as np
import matplotlib.pyplot as plt
import collections
#直方图均衡化
def histogram_equalization(histogram_e, lut_e, image_e):
sum_temp = 0
cf = []
for i in histogram_e:
sum_temp += i
cf.append(sum_temp)
for i, v in enumerate(lut_e):
lut_e[i] = int(255.0 * (cf[i] / sum_temp) + 0.5)
equalization_result = lut_e[image_e]
return equalization_result
# 计算灰度图的直方图
def draw_histogram(grayscale):
# 对图像进行通道拆分
hsi_i = grayscale[:, :, 2]
color_key = []
color_count = []
color_result = []
histogram_color = list(hsi_i.ravel()) # 将多维数组转换成一维数组
color = dict(collections.Counter(histogram_color)) # 统计图像中每个亮度级出现的次数
color = sorted(color.items(), key=lambda item: item[0]) # 根据亮度级大小排序
for element in color:
key = list(element)[0]
count = list(element)[1]
color_key.append(key)
color_count.append(count)
for i in range(0, 256):
if i in color_key:
num = color_key.index(i)
color_result.append(color_count[num])
else:
color_result.append(0)
color_result = np.array(color_result)
return color_result
x = []
for i in range(0, 256): # 横坐标
x.append(i)
# 原图及其直方图
origin = cv2.imread(r'D:\deng\smalljob\ImageSet\histeqColor.jpg')
original = cv2.cvtColor(origin,cv2.COLOR_BGR2RGB)
histogram_original = draw_histogram(original)
plt.bar(x, histogram_original) # 绘制原图直方图
plt.savefig('D:/deng/before_histogram.png')# 这里暂时不关闭绘制直方图的窗口,在处理图片后,再在该窗口绘制直方图作对比
before_histogram=cv2.imread(r'D:\deng\before_histogram.png')
lut = np.zeros(256, dtype=original.dtype) # 创建空的查找表
rgb_histogram = histogram_equalization(histogram_original, lut, original) # 均衡化处理
histogram_rgb_equalization = draw_histogram(rgb_histogram)
plt.bar(x, histogram_rgb_equalization) # 绘制原图直方图
plt.savefig('D:/deng/rgb_after_histogram.png')
plt.close()#关闭绘制直方图的窗口
rgb_after_histogram=cv2.imread(r'D:\deng\rgb_after_histogram.png')
#展示结果
plt.subplot(221),plt.imshow(original,cmap='gray')
plt.title('original')
plt.axis('off')
plt.subplot(222),plt.imshow(before_histogram,cmap='gray')
plt.title('before_histogram')
plt.axis('off')
plt.subplot(223),plt.imshow(rgb_histogram,cmap='gray')
plt.title('rgb_histogram')
plt.axis('off')
plt.subplot(224),plt.imshow(rgb_after_histogram,cmap='gray')
plt.title('rgb_after_histogram')
plt.axis('off')
plt.show()
运行结果:
以下是在HSI颜色模式下进行均衡化处理,即先将RGB图像转为HSI图像,再在此基础上进行均衡化,最后将图像转回RGB图像
#导入库
import cv2
import numpy as np
import matplotlib.pyplot as plt
import collections
from numpy import cos, arccos, sqrt, power, pi
def histogram_equalization(histogram_e, lut_e, image_e):
sum_temp = 0
cf = []
for i in histogram_e:
sum_temp += i
cf.append(sum_temp)
for i, v in enumerate(lut_e):
lut_e[i] = int(255.0 * (cf[i] / sum_temp) + 0.5)
equalization_result = lut_e[image_e]
return equalization_result
# 计算灰度图的直方图
def draw_histogram(grayscale):
# 对图像进行通道拆分
hsi_i = grayscale[:, :, 2]
color_key = []
color_count = []
color_result = []
histogram_color = list(hsi_i.ravel()) # 将多维数组转换成一维数组
color = dict(collections.Counter(histogram_color)) # 统计图像中每个亮度级出现的次数
color = sorted(color.items(), key=lambda item: item[0]) # 根据亮度级大小排序
for element in color:
key = list(element)[0]
count = list(element)[1]
color_key.append(key)
color_count.append(count)
for i in range(0, 256):
if i in color_key:
num = color_key.index(i)
color_result.append(color_count[num])
else:
color_result.append(0)
color_result = np.array(color_result)
return color_result
# HSI转RGB
def hsi_rgb(hsi):
if hsi.dtype.type == np.uint8:
hsi = (hsi).astype('float64') / 255.0
for k in range(hsi.shape[0]):
for j in range(hsi.shape[1]):
h, s, i = hsi[k, j, 0], hsi[k, j, 1], hsi[k, j, 2]
r, g, b = 0, 0, 0
if 0 <= h < 2/3*pi:
b = i * (1 - s)
r = i * (1 + s * cos(h) / cos(pi/3-h))
g = 3 * i - (b + r)
elif 2/3*pi <= h < 4/3*pi:
r = i * (1 - s)
g = i * (1 + s * cos(h-2/3*pi) / cos(pi - h))
b = 3 * i - (r + g)
elif 4/3*pi <= h <= 5/3*pi:
g = i * (1 - s)
b = i * (1 + s * cos(h - 4/3*pi) / cos(5/3*pi - h))
r = 3 * i - (g + b)
hsi[k, j, 0], hsi[k, j, 1], hsi[k, j, 2] = r, g, b
return (hsi * 255).astype('uint8')
# RGB转HSI
def rgb_hsi(rgb):
# 如果没有归一化处理,则需要进行归一化处理(传入的是[0,255]范围值)
if rgb.dtype.type == np.uint8:
rgb = rgb.astype('float64')/255.0
for i in range(rgb.shape[0]):
for j in range(rgb.shape[1]):
r, g, b = rgb[i, j, 0], rgb[i, j, 1], rgb[i, j, 2]
# 计算h
num = 0.5 * ((r-g)+(r-b))
den = sqrt(power(r-g, 2)+(r-b)*(g-b))
theta = arccos(num/den) if den != 0 else 0
rgb[i, j, 0] = theta if b <= g else (2*pi-theta)
# 计算s
rgb[i, j, 1] = (1 - 3 * min([r, g, b]) / (r+g+b)) if r+g+b != 0 else 0
# 计算i
rgb[i, j, 2] = 1 / 3 * (r+g+b)
return (rgb * 255).astype('uint8')
x = []
for i in range(0, 256): # 横坐标
x.append(i)
# 原图及其直方图
origin = cv2.imread(r'D:\deng\smalljob\ImageSet\histeqColor.jpg')
original = cv2.cvtColor(origin,cv2.COLOR_BGR2RGB)
histogram_original = draw_histogram(original)
plt.bar(x, histogram_original) # 绘制原图直方图
plt.savefig('D:/deng/before_histogram.png')
before_histogram=cv2.imread(r'D:\deng\before_histogram.png')
# rgb转hsi
#hsi_original = cv2.cvtColor(original, cv2.COLOR_RGB2HSV_FULL)
hsi_original = rgb_hsi(original)
# hsi在亮度分量上均衡化
histogram_hsi_original = draw_histogram(hsi_original)
plt.bar(x, histogram_hsi_original) # 绘制直方图
plt.savefig('D:/deng/hsi_before_histogram.png')
hsi_before_histogram=cv2.imread(r'D:\deng\hsi_before_histogram.png')
lut = np.zeros(256, dtype=hsi_original.dtype) # 创建空的查找表
hsi_histogram = histogram_equalization(histogram_hsi_original, lut, hsi_original) # 均衡化处理
# hsi转rgb
#result=cv2.cvtColor(hsi_histogram,cv2.COLOR_HSV2RGB_FULL)
result = hsi_rgb(hsi_histogram)
histogram_hsi_equalization = draw_histogram(result)
plt.bar(x, histogram_hsi_equalization) # 绘制原图直方图
plt.savefig('D:/deng/hsi_after_histogram.png')
plt.close()#关闭绘制直方图的窗口
hsi_after_histogram=cv2.imread(r'D:\deng\hsi_after_histogram.png')
#展示结果
plt.subplot(321),plt.imshow(original,cmap='gray')
plt.title('original')
plt.axis('off')
plt.subplot(322),plt.imshow(before_histogram,cmap='gray')
plt.title('before_histogram')
plt.axis('off')
plt.subplot(323),plt.imshow(hsi_original,cmap='gray')
plt.title('hsi_original')
plt.axis('off')
plt.subplot(324),plt.imshow(hsi_before_histogram,cmap='gray')
plt.title('hsi_before_histogram')
plt.axis('off')
plt.subplot(325),plt.imshow(result,cmap='gray')
plt.title('hsi_histogram')
plt.axis('off')
plt.subplot(326),plt.imshow(hsi_after_histogram,cmap='gray')
plt.title('hsi_after_histogram')
plt.axis('off')
plt.show()
运行结果:(图像出现绿色是因为在转换过程中出现了信息失真)