- pythonで「画像(行列)の整列」を実装したい人に向けた記事です。
- 相互相関係数を計算してシフト量を計算したのち、平行移動させて画像を整列させます。
- コードだけ載せます。
コード
#%%
import matplotlib.pyplot as plt
import numpy as np
import scipy
#%%
tx, ty = 0, 3
sizey, sizex = 11,15
#%%
xx,yy = np.meshgrid(np.arange(sizex),np.arange(sizey))
im0 = 1/(1+np.sqrt((xx-sizex/2)**2+(yy-sizey/2)**2))
im1 = 1/(1+np.sqrt(((xx-tx)-sizex/2)**2+((yy-ty)-sizey/2)**2))
fig,ax=plt.subplots(1,2)
ax[0].imshow(im0)
ax[1].imshow(im1)
ax[0].set_title("Original")
ax[1].set_title(f"Shifted ({tx},{ty})")
plt.plot()
#%%
# 偏差ベクトル
im0_dev = im0 -np.mean(im0)
im1_dev = im1 -np.mean(im1)
xcorr = scipy.signal.correlate2d(im0_dev,im1_dev, mode='same', boundary='fill', fillvalue=0)
# 相関係数に直す場合
xcorr /= (np.linalg.norm(im0_dev, ord=2) * np.linalg.norm(im1_dev, ord=2))
# 2Dの最大値の座標
y,x = np.unravel_index(np.argmax(xcorr), xcorr.shape)
plt.imshow(xcorr)
plt.title("Cross-correlation")
plt.plot(x,y, "x",c="r")
plt.colorbar()
plt.plot()
#%%
# 画像のシフト量
sy, sx = -(y-sizey//2), -(x-sizex//2)
# 平行移動
mat = [[1,0,sy],
[0,1,sx],
[0,0,1]]
im2 = scipy.ndimage.interpolation.affine_transform(im1, mat, order=1)
fig,ax=plt.subplots(1,3, sharey=True, tight_layout=True)
ax[0].imshow(im0)
ax[1].imshow(im1)
ax[2].imshow(im2,vmax=im0.max(), vmin=im0.min())
ax[0].set_title("Original")
ax[1].set_title(f"Shifted ({tx},{ty})")
ax[2].set_title("Aligned")
plt.plot()※ #%%はVScodeのセル区切りです。
使い方
シフトさせたダミー画像の用意
ダミー画像として、y方向にシフトさせた画像を用意します。ｙ方向に３動かしました。
tx, ty = 0, 3
sizey, sizex = 11,15
#%%
xx,yy = np.meshgrid(np.arange(sizex),np.arange(sizey))
im0 = 1/(1+np.sqrt((xx-sizex/2)**2+(yy-sizey/2)**2))
im1 = 1/(1+np.sqrt(((xx-tx)-sizex/2)**2+((yy-ty)-sizey/2)**2))
これらの画像の二次元の相互相関係数を計算することで、どれだけシフトさせたかを計算します。
二次元相互相関によるシフト量計算
二次元の相互相関はnumpyには用意されていません。scipyのscipy.signal.correlate2dを用います。
この関数はcorrelateと言いつつも何も加工をしてくれないので、手動で偏差ベクトル（平均値を減算）にして前処理します。
また、今回は必須ではありませんが、相互相関係数（0から１までの係数）にする場合はやはり手動で直す必要があります。例のようにベクトルの大きさで割ってください。
二次元相互相関係数xcorrが求まったら、最大値を探します。二次元だとnp.argmax()だけでは座標が得られないので、面倒ですがnp.unravel_index()と組み合わせます。
# 偏差ベクトル
im0_dev = im0 -np.mean(im0)
im1_dev = im1 -np.mean(im1)
xcorr = scipy.signal.correlate2d(im0_dev,im1_dev, mode='same', boundary='fill', fillvalue=0)
# 相関係数に直す場合
xcorr /= (np.linalg.norm(im0_dev, ord=2) * np.linalg.norm(im1_dev, ord=2))
# 2Dの最大値の座標
y,x = np.unravel_index(np.argmax(xcorr), xcorr.shape)
アフィン変換による整列
以上の計算したシフト量に基づいて、画像をシフトさせます。
画像のシフトは色々ありますが、ここではscipy.ndimage.interpolation.affine_transform()を用いました。似たような関数がopencv（cv2）にもありますが、scipyで統一するのでこちらを使います。
スライスで切り取るやり方と違って、はみ出てエラーになる心配がなくとても便利です。
引数のmatにはアフィン変換の行列(同次変換の行列)を指定しますが、ここはシフトだけなので簡単ですね。
また、ここではシフト量が整数なので意味はないのですが、order=1は補完方法を線形(一次)にしています。
シフト量に少数含む場合ときにサブピクセル精度で並行シフトするのに使えます。
Alignedの画像を見れば、ShiftedがOriginalに整列されたことが確認できました。
# 画像のシフト量
sy, sx = -(y-sizey//2), -(x-sizex//2)
# 平行移動
mat = [[1,0,sy],
[0,1,sx],
[0,0,1]]
im2 = scipy.ndimage.interpolation.affine_transform(im1, mat, order=1)
コメント