直方图均衡化
-
灰度映射函数：
s = T ( r ) , 0 ≤ r ≤ L − 1 (3.8) s = T(r), \quad 0\leq r \leq L -1 \tag{3.8} s=T(r),0≤r≤L−1(3.8) -
假设：
-
(1) T ( r ) T(r) T(r)在区间 0 ≤ r ≤ L − 1 0 \leq{r} \leq{L-1} 0≤r≤L−1 上是一个单调递增函数。
-
(2) 对于 0 ≤ r ≤ L − 1 0 \leq{r} \leq{L-1} 0≤r≤L−1 ，有 0 ≤ T ( r ) ≤ L − 1 0 \leq{T(r)} \leq{L-1} 0≤T(r)≤L−1 。
-
逆变换
r = T − 1 ( s ) , 0 ≤ s ≤ L − 1 (3.9) r = T^{-1}(s), \quad 0 \leq s \leq L-1 \tag{3.9} r=T−1(s),0≤s≤L−1(3.9)
将条件(1)改为：
- (1’) T ( r ) T(r) T(r)在区间 0 ≤ r ≤ L − 1 0 \leq{r} \leq{L-1} 0≤r≤L−1 上是一个严格单调递增函数。
(1)中
T
(
r
)
T(r)
T(r)单调递增的条件保证输出灰度值从不小于对应的输入值，从而防止灰度反转而产生伪像。
(2)保证输出灰度的范围与输入的范围相同。
(1’)保证从
s
s
s返回到
r
r
r的映射是一对一的，从而防止出现歧义。
图像的灰度可视为区间
[
0
,
L
−
1
]
[0, L-1]
[0,L−1]内的一个随机变量。令
p
r
(
r
)
p_{r}(r)
pr(r)和
p
s
(
s
)
p_{s}(s)
ps(s)表示两幅不同图像中灰度值
r
r
r和
s
s
s的PDF（概率密度函数）。
p
p
p的下标表明
p
r
(
r
)
p_{r}(r)
pr(r)和
p
s
(
s
)
p_{s}(s)
ps(s)是不同的函数。若已知
p
r
(
r
)
p_{r}(r)
pr(r)和
T
(
r
)
T(r)
T(r)，且
T
(
r
)
T(r)
T(r)是连续的且在感兴趣的值域上是可微的，则变换（映射）后的变量
s
s
s的PDF是
p
s
(
s
)
=
p
r
(
r
)
∣
d
r
d
s
∣
(3.10)
p_{s}(s) = p_{r}(r) \bigg\lvert \frac{\text{d}r}{\text{d}s}\bigg\rvert \tag{3.10}
ps(s)=pr(r)∣∣∣∣dsdr∣∣∣∣(3.10)
可以看到输出灰度变量
s
s
s的PDF是由输入灰度的PDF和所有的变换函数确定的
图像处理中的一个特别重要的变换函数是
s
=
T
(
r
)
=
(
L
−
1
)
∫
0
r
p
r
(
w
)
d
w
(3.11)
s = T(r) = (L-1) \int_{0}^{r} p_r(w) \text{d} w \tag{3.11}
s=T(r)=(L−1)∫0rpr(w)dw(3.11)
w
w
w是一个假积分变量，右侧的积分是随机变量
r
r
r的累积分布函数（CDF）
根据莱布尼茨积分法则可知
d
r
d
s
=
d
T
(
r
)
d
r
=
(
L
−
1
)
d
d
r
[
∫
0
r
p
r
(
w
)
d
w
]
=
(
L
−
1
)
p
r
(
r
)
(3.12)
\frac{\text{d}r}{\text{d}s} = \frac{\text{d}T(r)}{\text{d}r} = (L - 1) \frac{\text{d}}{\text{d}r} \Bigg[\int_{0}^{r} p_r(w) \text{d} w \Bigg] = (L-1)p_{r}(r)\tag{3.12}
dsdr=drdT(r)=(L−1)drd[∫0rpr(w)dw]=(L−1)pr(r)(3.12)
用这个结果代替式（3.10）的
d
r
d
s
\frac{\text{d}r}{\text{d}s}
dsdr，并且所有的概率值都是正的，有：
p
s
(
s
)
=
p
r
(
r
)
∣
d
r
d
s
∣
=
p
r
(
r
)
∣
1
(
L
−
1
)
p
r
(
r
)
∣
=
1
L
−
1
,
0
≤
s
≤
L
−
1
(3.10)
p_{s}(s) = p_{r}(r) \bigg\lvert \frac{\text{d}r}{\text{d}s}\bigg\rvert = p_{r}(r) \bigg\lvert \frac{1}{(L-1)p_{r}(r)}\bigg\rvert = \frac{1}{L - 1}, \;0 \leq s \leq L - 1\tag{3.10}
ps(s)=pr(r)∣∣∣∣dsdr∣∣∣∣=pr(r)∣∣∣∣(L−1)pr(r)1∣∣∣∣=L−11,0≤s≤L−1(3.10)
直方图均衡化或直方图线性变换
式（3.11)变换的离散形式为：
s
k
=
T
(
r
k
)
=
(
L
−
1
)
∑
j
=
0
k
p
r
(
r
j
)
,
k
=
0
,
1
,
2
,
…
,
L
−
1
(3.15)
s_{k} = T(r_{k}) = (L -1) \sum_{j=0}^k p_{r}(r_{j}),\quad k = 0, 1, 2, \dots, L-1 \tag{3.15}
sk=T(rk)=(L−1)j=0∑kpr(rj),k=0,1,2,…,L−1(3.15)
def my_calhist(img):
"""
histogram equalization
param: input img: uint8[0, 255] grayscale image
return uint8[0, 255] grayscale image after histogram equalization
"""
hist, bins = my_hist(img, bins=256, normalized=True)
#--------------------------------Numpy-------------------
hist_cumsum = np.round(np.cumsum(hist) * 255).astype(int)
img_dst = img.copy()
img_dst = hist_cumsum[img_dst]
#--------------------------------loop----------------------
# map_dict = {}
# map_dict = map_dict.fromkeys(bins.astype(int), 0)
# for i in bins:
# s = np.round(255 * hist[:i].sum()).astype(int)
# map_dict[i] = s
# height, width = img.shape[:2]
# img_dst = np.zeros([height, width], np.uint8)
# for h in range(height):
# for w in range(width):
# img_dst[h, w] = map_dict[img[h, w]] # dict 用[ ]
# img_dst = np.clip(img_dst, 0, 255).astype(np.uint8)
# return img_dst, map_dict
return img_dst, _
# 直方图均衡化
import time
img_1st = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0316(1)(top_left).tif', 0)
img_2nd = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0316(2)(2nd_from_top).tif', 0)
img_3rd = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0316(3)(third_from_top).tif', 0)
img_4th = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0316(4)(bottom_left).tif', 0)
img_list =['img_1st', 'img_2nd', 'img_3rd', 'img_4th']
fig = plt.figure(figsize=(20, 20))
start_time = time.time()
for i in range(len(img_list)):
# plot orginal image
ax = fig.add_subplot(4, 4, i+1)
ax.imshow(eval(img_list[i]), cmap='gray', vmin=0, vmax=255), ax.set_title('Original')
# plot histogram equalization image
img_dst, _ = my_calhist(eval(img_list[i]))
ax1 = fig.add_subplot(4, 4, i+5)
ax1.imshow(img_dst, cmap='gray', vmin=0, vmax=255), ax1.set_title('After Histogram Equalization')
# plot orginal hist
ax2 = fig.add_subplot(4, 4, i+9)
hist, bins = my_hist(eval(img_list[i]), bins=256)
ax2.bar(bins, hist), ax2.set_title('Original Histogram')
# plot hist after histogram equalization
ax3 = fig.add_subplot(4, 4, i+13)
hist, bins = my_hist(img_dst, bins=256)
ax3.bar(bins, hist), ax3.set_title('Hist After Histogram Equalization')
elapse = time.time() - start_time
print(f'Elapse -> {elapse}s')
plt.tight_layout()
plt.show()
Elapse -> 2.3357810974121094s
def my_calhist(img):
"""
histogram equalization
param: input img: uint8[0, 255] grayscale image
return uint8[0, 255] grayscale image after histogram equalization
"""
hist, bins = my_hist(img, bins=256, normalized=True)
hist_cumsum = np.round(np.cumsum(hist) * 255).astype(int)
height, width = img.shape[:2]
img_dst = np.zeros([height, width], np.uint8)
for h in range(height):
for w in range(width):
img_dst[h, w] = hist_cumsum[img[h, w]] # dict 用[ ]
img_dst = np.clip(img_dst, 0, 255).astype(np.uint8)
return img_dst, hist_cumsum
# 直方图均衡化转换函数
img_1st = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0316(1)(top_left).tif', 0)
img_2nd = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0316(2)(2nd_from_top).tif', 0)
img_3rd = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0316(3)(third_from_top).tif', 0)
img_4th = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0316(4)(bottom_left).tif', 0)
img_list =['img_1st', 'img_2nd', 'img_3rd', 'img_4th']
fig = plt.figure(figsize=(6, 6))
for i in range(len(img_list)):
# plot histogram equalization transform
_, map_dict = my_calhist(eval(img_list[i]))
ax1 = fig.gca()
ax1.plot(map_dict)
plt.xlabel('Original gray value', fontsize=12)
plt.ylabel('Transformed gray value', fontsize=12)
plt.xticks([0, 64, 128, 192, 255])
plt.yticks([0, 64, 128, 192, 255])
plt.xlim([0, 255])
plt.ylim([0, 255])
plt.tight_layout()
plt.show()
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
