互信息与相关性的图像配准
一 基于互信息与相关性的图像配准
1 互信息定义及图像配准
互信息定义如下：
\[MI(X;Y) = \sum\limits_{x}\sum\limits_{y}p_{(X,Y)}(x,y)\log\cfrac{p_{(X,Y)}(x,y)}{{p_X(x)}\cdot{p_Y(y)}}
\]
从定义中可以看出，当X和Y两个变量趋近于独立时，
\[p_{(X,Y)}(x,y)={p_X(x)}\cdot{p_Y(y)}
\]
此时，
\[\log\cfrac{p_{(X,Y)}(x,y)}{{p_X(x)}\cdot{p_Y(y)}}\approx 0，MI(X;Y)有最小值。
\]
反之，若X和Y两个变量越相关，
\[MI(X;Y)则越大。
\]
基于以上理论，设有模板图片fixed、待配准图片moving，若
\[fixed ( x , y ) = moving( x − △ x , y − △ y )
\]
则当moving平移 △x、 △y 时，以fixed、moving图像二维亮度直方图(表示联合分布)作为输入，求得的互信息MI有最大值。
2 相关性的定义及图像配准
皮尔逊相关系数：
\[ ρX,Y=\cfrac{cov(X,Y)}{σXσY}=\cfrac{[(X−μX)(Y−μY)]}{σXσYE}
\]
\[ρX,Y =\cfrac{\sum\limits_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{{\sum\limits_{i=1}^n(x_i-\bar{x}} )^2}\sqrt{{\sum\limits_{i=1}^n(y_i-\bar{y}} )^2}}\\ \\ \quad\quad\ = \cfrac{n\sum x_iy_i-\sum x_i \sum y_i}{\sqrt{n \sum x_i^2-(\sum x_i)^2}\sqrt{n \sum y_i^2-(\sum y_i)^2}}
\]
基于以上理论，设有模板图片fixed、待配准图片moving，若
\[fixed ( x , y ) = moving( x − △ x , y − △ y )
\]
则当moving平移 △x、 △y 时，将fixed、moving图像按照一维数据排列，进行相关性计算，求得的相关性有最大值；相关性越接近1，配准效果越好，反之，相关性越接近0，配准效果越差。
二 模板图像fixed、（待）配准图像moving的可视化
1 fixed 和 moving 图像的直方图
注：配准后图像的缺失边界可以用fixed图像对应的像素来填充，这样差分图对应的像素值就为0。
2 fixed 和 moving 图像的散点图
从散点图上可以看出，fixed 和 待配准的 moving 图像的像素分布是相互独立的，而fixed 和 配准的 moving 图像的像素分布表现出一定的线性关系。
3 fixed 和 moving 图像的2D直方图
同样，从2D直方图上可以看出，fixed 和 待配准的 moving 图像的像素分布线性度不好，而fixed 和 配准的 moving 图像的像素分布表现出一定的线性关系。。
4 fixed 和 moving 图像的差分图
三 图像配准的代码
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]))
def correlation_information(img_fixed, img_moved):
def NCC(img1, img2):
# # 图片2的标准差
# print(np.std(img2))
# 相关系数，这里使用的是有偏估计
# 协方差 Cov(X, Y) = Σ((Xᵢ - μₓ)(Yᵢ - μᵧ)) / n
# 方差 D(X) = (∑(x - μ)²) / n
# R = Con(X, Y) / (sqrt(D(X)) * sqrt(D(Y)))
con = np.mean(np.multiply((img1 - np.mean(img1)), (img2 - np.mean(img2)))) # 协方差
img1Std = np.std(img1)
img2Std = np.std(img2)
NCCValue = con / (img1Std * img2Std)
return NCCValue
r, _ = stats.pearsonr(img_fixed.ravel(), img_moved.ravel())
# ncc_cal_in = np.corrcoef(img_fixed.astype(np.float32).flatten(), img_moved.astype(np.float32).flatten())[0][1]
# ncc_cal = NCC(img_fixed, img_moved)
# print(r, ncc_cal_in, ncc_cal)
return r
METHOD_REGISTRATION_MI = 0
METHOD_REGISTRATION_CI = 1
def ImageRegistration(img_fixed, img_moved, max_offset_XY=(5, 5), method=METHOD_REGISTRATION_MI):
large = -np.inf
large_ij = []
h, w = img_moved.shape[:2]
ox, oy = max_offset_XY
dirs = [(-1, -1), (-1, 1), (1, 1), (1, -1)]
for i in range(oy):
for j in range(ox):
if i == 0 and j == 0:
ret = 0
if method == METHOD_REGISTRATION_MI:
hist_2d_moved, x_edges, y_edges = np.histogram2d(
img_fixed.ravel(),
img_moved.ravel(),
bins=255, range=[[0, 255], [0, 255]])
ret = mutual_information(hist_2d_moved)
if method == METHOD_REGISTRATION_CI:
ret = correlation_information(img_fixed, img_moved)
if ret > large:
large = ret
large_ij = [i, j]
continue
for _dir in dirs:
deltaY = i * _dir[1]
deltaX = j * _dir[0]
M = np.array([[1, 0, deltaX], [0, 1, deltaY]], dtype=np.float)
img_moved_modify = cv2.warpAffine(img_moved, M, dsize=(w, h))
ret = 0
if method == METHOD_REGISTRATION_MI:
hist_2d_moved, x_edges, y_edges = np.histogram2d(
img_fixed.ravel(),
img_moved_modify.ravel(),
bins=256, range=[[0, 255], [0, 255]])
ret = mutual_information(hist_2d_moved)
if method == METHOD_REGISTRATION_CI:
ret = correlation_information(img_fixed, img_moved_modify)
if ret > large:
large = ret
large_ij = [deltaY, deltaX]
offset_y, offset_x = large_ij
print("ret = {}, offset_y = {}, offset_x = {} pixel".format(large, offset_y, offset_x))
M = np.array([[1, 0, offset_x], [0, 1, offset_y]], dtype=np.float)
img_moved_final = cv2.warpAffine(img_moved, M, dsize=(w, h))
return img_moved_final, offset_x, offset_y
