直方图匹配（规定化）

• 连续灰度
  $$s=T(r)=(L-1){\int }_0^r{p}_r(w)\text{d}w\text{\hspace{1em}(3.17)}$$
  定义关于变量$z$的一个函数$G$，它具有如下性质：
  $$G(z)=(L-1){\int }_0^z{p}_z(v)\text{d}v=s\text{\hspace{1em}(3.18)}$$
  因此$z$必定满足条件：
  $$z=G^{-1}(s)=G^{-1}[T(r)]\text{\hspace{1em}(3.19)}$$

使用如下步骤可以得到一幅灰度级具有规定PDF的图像：

1. 由输入图像得到式（3.17）中使用的$p_r(r)$。
2. 在式（3.18）中使用规定的PDF即$p_z(z)$得到函数$G(z)$。
3. 计算反变换$z=G^{-1}(s)$；输出图像中的像素值是$s$。对于均衡化后的图像中值为$s$的每个像素执行逆映射$z=G^{-1}(s)$，得到输出图像中的对应像素。使用这个变换处理完所有像素后，输出图像的PDF即$p_z(z)$将等于规定的PDF。

• 离散灰度
  $$s_k=T(r_k)=(L-1)\sum _{j=0}^k{p}_r(r_j),k=0,1,2,\dots ,L-1\text{\hspace{1em}(3.20)}$$
  $$G(z_q)=(L-1)\sum _{i=0}^q{p}_z(z_i),q=0,1,2,\dots ,L-1\text{\hspace{1em}(3.21)}$$
  $$G(z_q)=s_k\text{\hspace{1em}(3.22)}$$
  $$z_q=G^{-1}(s_k)\text{\hspace{1em}(3.23)}$$

离散直方图规定化的过程

1. 计算输入图像的直方图$p_r(r)$，并在式（3.20）中用它将输入图像中的灰度映射到直方图均衡化后的图像中的灰度。将得到的值$s_k$四舍五入到整数区间$[0,L-1]$。
2. 用式（3.21）对$q=0,1,2,\dots ,L-1$计算函数$G(z_q)$的所有值，其中$p_z(z_i)$是规定直方图的值。将$G$的值四舍五入到区间$[0,L-1]$内的整数。并存储到一个查找表中。
3. 对$s_k,k=0,1,2,\dots ,L-1$的每个值，用步骤2中存储的$G$值找到$z_q$的对应值，使得$G(z_q)$最接近$s_k$。存储从$s$到$z$的这些映射。当多个$z_q$值给出相同的匹配（即映射不唯一）时，按照约定选择最小的值。
4. 使用步骤3中找到的映射，将每个均衡化后的像素$s_k$值映射到直方图规定化图像中值为$z_q$的对应像素，并形成直方图规定化后的图像。

def my_histogram_matching(img_src, img_dst):"""historgram specification for two input imagesparam: input img_src: input image used for histogram specification , uint8[0, 255] grayscale imageparam: input img_dst: input image for histogram specification , uint8[0, 255] grayscale imagereturn: uint8[0, 255] grayscale image after histogram specification """bins = 256# img_src -> input, img_dst -> matchinghist_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 = 1temp_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()