图像配准之二：基于互信息_图像配准互信息

1.互信息定义

互信息定义如下：

$MI(X;Y)=\sum _x\sum _y{p}_{(X,Y)}(x,y)\log \frac{p_{(X,Y)}(x,y)}{p_X(x)\cdot p_Y(y)}$

从定义中可以看出，当X和Y两个变量趋近于独立时，$p_{(X,Y)}(x,y)\approx p_X(x)\cdot p_Y(y)$。

此时，$\log \frac{p_{(X,Y)}(x,y)}{p_X(x)\cdot p_Y(y)}\approx 0$，$MI(X;Y)$有最小值。

反之，若X和Y两个变量越相关，$MI(X;Y)$则越大。

更多内容可参考KL散度 JS散度 熵 互信息。

2.基于互信息的图像配准

基于以上理论，设有模板图片$I_1$、待配准图片$I_2$，若$I_1(x,y)=I_2(x-△x,y-△y)$，

则当$I_2$平移$x,y$ 时，以$I_1$、$I_2$图像二维亮度直方图(表示联合分布)作为输入，求得的互信息MI有最大值。

2.1 模板图片$I_1$、待配准图片$I_2$的可视化

同图像配准之一：FFT，这里也选择lena图片作为输入。

• 模板图片$I_1$、待配准图片$I_2$各自的直方图信息

截取两块区域作为输入，

img = plt.imread("lena.bmp")
img1=img[10:,20:]
img2=img[:-10,:-20]
fig, axes = plt.subplots(2, 2)
axes[0,0].imshow(img1)
axes[0,1].imshow(img2)

axes[1,0].hist(img1.ravel(), bins=20)
axes[1,0].set_title('img1 histogram')
axes[1,1].hist(img2.ravel(), bins=20)
axes[1,1].set_title('img2 histogram')

plt.show(block=True)
1
2
3
4
5
6
7
8
9
10
11
12
13

• 模板图片与待配准图片像素的散点图

plt.plot(img1.ravel(), img2.ravel(), '.')
plt.xlabel('img1 signal')
plt.ylabel('img2 signal')
plt.title('img1 vs img2 signal')
plt.show(block=True)
1
2
3
4
5

从这个图上也可以看出，两张图像素分布相对独立。

• 模板图片与待配准图片的2D直方图

hist_2d, x_edges, y_edges = np.histogram2d(img1.ravel(),img2.ravel(),bins=20)
plt.imshow(hist_2d.T, origin='lower')
plt.xlabel('img1 signal bin')
plt.ylabel('img2 signal bin')
plt.show(block=True)
1
2
3
4
5

2.2 将2D直方图作为输入，计算互信息

2.2.1 通过2D直方图计算互信息

def mutual_information(hgram):
    """ Mutual information for joint histogram
    """
    # Convert bins counts to probability values
    pxy = hgram / float(np.sum(hgram))
    px = np.sum(pxy, axis=1) # marginal for x over y
    py = np.sum(pxy, axis=0) # marginal for y over x

    px_py = px[:, None] * py[None, :] # Broadcast to multiply marginals
    nzs = pxy > 0 # Only non-zero pxy values contribute to the sum

    return np.sum(pxy[nzs] * np.log2(pxy[nzs] / px_py[nzs]))
1
2
3
4
5
6
7
8
9
10
11
12

2.2.2 遍历两个方向的偏移，计算互信息最大值

large=-np.inf
larg_ij=[]
w,h=img2.shape


for i in range(30):
    for j in range(30):

        img2_moved = np.zeros(img2.shape)
        img2_moved[:w-i, :h-j] = img2[i:, j:]

        hist_2d_moved, x_edges, y_edges = np.histogram2d(
            img1.ravel(),
            img2_moved.ravel(),
            bins=20)

        mi=mutual_information(hist_2d_moved)
        if mi>large:
            large=mi
            larg_ij=[i,j]
print(large,larg_ij)
#2.793162258691581 [10, 20]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

得到正确的偏移量 10, 20。

2.3 根据偏移量，再次可视化1D、2D直方图

• 模板图片、配准后图片的1D直方图

fig, axes = plt.subplots(2, 2)
img2_moved = np.zeros(img2.shape)
i,j=10,20
img2_moved[:w-i, :h-j] = img2[i:, j:]
axes[0,0].imshow(img1)
axes[0,1].imshow(img2_moved)

axes[1,0].hist(img1.ravel(), bins=20)
axes[1,0].set_title('img1 histogram')
axes[1,1].hist(img2_moved.ravel(), bins=20)
axes[1,1].set_title('img2_moved histogram')

plt.show(block=True)

1
2
3
4
5
6
7
8
9
10
11
12
13
14

• 模板图片与配准后图片的像素散点图

plt.plot(img1.ravel(), img2_moved.ravel(), '.')
plt.xlabel('img1 signal')
plt.ylabel('img2_moved signal')
plt.title('img1 vs img2_moved signal')
1
2
3
4

可见，两张图片像素分布相关性很强。

• 模板图片与配准后图片的2D直方图

hist_2d, x_edges, y_edges = np.histogram2d(img1.ravel(),img2_moved.ravel(),bins=20)
plt.imshow(hist_2d.T, origin='lower')
plt.xlabel('img1 signal bin')
plt.ylabel('img2_moved signal bin')
plt.show(block=True)
1
2
3
4
5

参考文献

[1] Mutual information as an image matching metric
[2] wiki/Mutual_information