目錄
一 基於互資訊與相關性的影像配準
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