第3章 Python 数字图像处理(DIP) - 灰度变换与空间滤波9 - 直方图处理 - 直方图匹配（规定化）灰度图像，彩色图像都适用

直方图匹配（规定化）

• 连续灰度
  $$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 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()