直方图匹配(规定化)
- 连续灰度
s = T ( r ) = ( L − 1 ) ∫ 0 r p r ( w ) d w (3.17) s = T(r) = (L-1) \int_{0}^{r} p_r(w) \text{d} w \tag{3.17} s=T(r)=(L−1)∫0rpr(w)dw(3.17)
定义关于变量 z z z的一个函数 G G G,它具有如下性质:
G ( z ) = ( L − 1 ) ∫ 0 z p z ( v ) d v = s (3.18) G(z) = (L - 1) \int_{0}^{z} p_z(v) \text{d} v = s \tag{3.18} G(z)=(L−1)∫0zpz(v)dv=s(3.18)
因此 z z z必定满足条件:
z = G − 1 ( s ) = G − 1 [ T ( r ) ] (3.19) z = G^{-1}(s) = G^{-1}[T(r)] \tag{3.19} z=G−1(s)=G−1[T(r)](3.19)
使用如下步骤可以得到一幅灰度级具有规定PDF的图像:
- 由输入图像得到式(3.17)中使用的 p r ( r ) p_{r}(r) pr(r)。
- 在式(3.18)中使用规定的PDF即 p z ( z ) p_{z}(z) pz(z)得到函数 G ( z ) G(z) G(z)。
- 计算反变换 z = G − 1 ( s ) z=G^{-1}(s) z=G−1(s);输出图像中的像素值是 s s s。对于均衡化后的图像中值为 s s s的每个像素执行逆映射 z = G − 1 ( s ) z=G^{-1}(s) z=G−1(s),得到输出图像中的对应像素。使用这个变换处理完所有像素后,输出图像的PDF即 p z ( z ) p_{z}(z) pz(z)将等于规定的PDF。
- 离散灰度
s k = T ( r k ) = ( L − 1 ) ∑ j = 0 k p r ( r j ) , k = 0 , 1 , 2 , … , L − 1 (3.20) 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.20} sk=T(rk)=(L−1)j=0∑kpr(rj),k=0,1,2,…,L−1(3.20)
G ( z q ) = ( L − 1 ) ∑ i = 0 q p z ( z i ) , q = 0 , 1 , 2 , … , L − 1 (3.21) G(z_{q}) = (L -1) \sum_{i=0}^q p_{z}(z_{i}),\quad q = 0, 1, 2, \dots, L-1 \tag{3.21} G(zq)=(L−1)i=0∑qpz(zi),q=0,1,2,…,L−1(3.21)
G ( z q ) = s k (3.22) G(z_{q}) = s_{k} \tag{3.22} G(zq)=sk(3.22)
z q = G − 1 ( s k ) (3.23) z_{q} = G^{-1}(s_{k}) \tag{3.23} zq=G−1(sk)(3.23)
离散直方图规定化的过程
- 计算输入图像的直方图 p r ( r ) p_{r}(r) pr(r),并在式(3.20)中用它将输入图像中的灰度映射到直方图均衡化后的图像中的灰度。将得到的值 s k s_{k} sk四舍五入到整数区间 [ 0 , L − 1 ] [0, L-1] [0,L−1]。
- 用式(3.21)对 q = 0 , 1 , 2 , … , L − 1 q = 0, 1, 2, \dots, L-1 q=0,1,2,…,L−1计算函数 G ( z q ) G(z_{q}) G(zq)的所有值,其中 p z ( z i ) p_{z}(z_{i}) pz(zi)是规定直方图的值。将 G G G的值四舍五入到区间 [ 0 , L − 1 ] [0, L-1] [0,L−1]内的整数。并存储到一个查找表中。
- 对 s k , k = 0 , 1 , 2 , … , L − 1 s_{k}, k=0, 1, 2, \dots, L-1 sk,k=0,1,2,…,L−1的每个值,用步骤2中存储的 G G G值找到 z q z_{q} zq的对应值,使得 G ( z q ) G(z_{q}) G(zq)最接近 s k s_{k} sk。存储从 s s s到 z z z的这些映射。当多个 z q z_{q} zq值给出相同的匹配(即映射不唯一)时,按照约定选择最小的值。
- 使用步骤3中找到的映射,将每个均衡化后的像素 s k s_{k} sk值映射到直方图规定化图像中值为 z q z_{q} zq的对应像素,并形成直方图规定化后的图像。
def my_histogram_matching(img_src, img_dst):
"""
historgram specification for two input images
param: input img_src: input image used for histogram specification , uint8[0, 255] grayscale image
param: input img_dst: input image for histogram specification , uint8[0, 255] grayscale image
return: uint8[0, 255] grayscale image after histogram specification
"""
bins = 256
# img_src -> input, img_dst -> matching
hist_src, bins_src = my_hist(img_src, bins=256, normalized=True)
hist_dst, bins_dst = my_hist(img_dst, bins=256, normalized=True)
# [原灰度级,均衡化的值]
list_src = list([x, y] for x in np.arange(bins) for y in np.arange(bins) if y == 0)
for i in bins_src:
s = np.round(255 * hist_src[:i].sum()).astype(int)
list_src[i][1] = s
# [原灰度级,均衡化的值]
list_dst = list([x, y] for x in np.arange(bins) for y in np.arange(bins) if y == 0)
for i in bins_dst:
s = np.round(255 * hist_dst[:i].sum()).astype(int)
list_dst[i][1] = s
# 映射的关系列表初始化,以下映射关系算法
list_dst_1 = list([x, y] for x in np.arange(bins) for y in np.arange(bins) if y == 0)
for i in range(256):
minv = 1
temp_dst = list_dst[i][1]
for j in range(256):
temp_src = list_src[j][1]
if abs(temp_dst - temp_src) < minv:
minv = abs(temp_dst - temp_src)
idx = int(j)
# print(f'dst -> {temp_dst}, src -> {temp_src}, idx -> {idx}')
list_dst_1[i][1] = idx
#------------------------------Numpy-------------
map_dict = np.array(list_dst_1)[:, 1]
img_result = img_dst.copy()
img_result = map_dict[img_result]
#------------------------------loop-----------------
# height, width = img_dst.shape[:2]
# img_result = np.zeros([height, width], np.uint8)
# for h in range(height):
# for w in range(width):
# img_result[h, w] = list_dst_1[img_dst[h, w]][1]
# img_result = np.clip(img_result, 0, 255).astype(np.uint8)
return img_result
# Grayscale image histogram matching, different shape image works
img_1 = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0307(a)(intensity_ramp).tif', 0)
img_2 = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH03/Fig0316(4)(bottom_left).tif', 0)
hist_1, bins_1 = my_hist(img_1, bins=256, normalized=True)
hist_2, bins_2 = my_hist(img_2, bins=256, normalized=True)
img_dst = my_histogram_matching(img_1, img_2,)
hist_dst, bins_dst = my_hist(img_dst, bins=256, normalized=True)
plt.figure(figsize=(15, 10))
plt.subplot(2, 3, 1), plt.imshow(img_1, 'gray', vmin=0, vmax=255), plt.title('src')
plt.subplot(2, 3, 2), plt.imshow(img_2, 'gray', vmin=0, vmax=255), plt.title('dst')
plt.subplot(2, 3, 3), plt.imshow(img_dst, 'gray', vmin=0, vmax=255), plt.title('Result')
plt.subplot(2, 3, 4), plt.bar(bins_1, hist_1), plt.xlim([0, 255])
plt.subplot(2, 3, 5), plt.bar(bins_2, hist_2), plt.xlim([0, 255])
plt.subplot(2, 3, 6), plt.bar(bins_dst, hist_dst), plt.xlim([0, 255])
plt.tight_layout()
plt.show()
# RGB 直方图匹配 opencv读入的是BGR图像,[..., ::-1]可以把BGR转为RGB的图像
img_1 = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH06/Fig0651(a)(flower_no_compression).tif',)[..., ::-1]
img_2 = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH06/Fig0638(a)(lenna_RGB).tif', )[..., ::-1]
img_dst = np.zeros_like(img_2, dtype=np.uint8)
for i in range(3):
img_dst[..., i] = my_histogram_matching(img_1[..., i], img_2[..., i])
plt.figure(figsize=(15, 10))
plt.subplot(2, 3, 1), plt.imshow(img_1, vmin=0, vmax=255), plt.title('src')
plt.subplot(2, 3, 2), plt.imshow(img_2, vmin=0, vmax=255), plt.title('dst')
plt.subplot(2, 3, 3), plt.imshow(img_dst, vmin=0, vmax=255), plt.title('Result')
plt.tight_layout()
plt.show()