Help us understand the problem. What is going on with this article?

画像重ね合わせプログラム(更新中)

1. はじめに

本稿では平行移動、回転、スケールされた画像(slave)を元の画像(master)に重ね合わせるpythonのプログラムを解説する。
Figure_1-1.png
Figure_2.png

まず、プログラムの動作環境について説明しておく。プログラムはLinux/Mac/Windows10のAnanconda/Python3.6で動作するので、Anacondaをインストールしていない場合、まずは以下からAnacondaをダウンロードしてきて、インストールする。
https://www.anaconda.com/download/

Windows10を前提に話を進める。command promptを立ち上げて次のコマンドを実行する。

$ conda create -n py36 python=3.6
Proceed ([y]/n)? y
$ activate py36
$ pip install numpy opencv-python matplotlib

(Macの場合はmatplotlibのインストールでもう一手間かかる)
https://qiita.com/itoru257/items/c52cab383eaae4a242d4

これで、今回のプログラムを動かすpython3.6の環境が構築できた。

プログラムソースコードをダウロードして解凍する。
https://github.com/logicool-repo/fft-based_image_coregistration/archive/master.zip

さきほどのコマンドプロンプトでcdを解凍プログラムのディレクトリに移動させ、以下のコマンドを叩いてプログラムを実行する。lenaのサンプル画像の重ね合わせ結果が出力されるはず。

$ python main.py

2. 重ね合わせアルゴリズムの紹介

[1] B.S. Reddy and B.N. Chatterji, "An FFT-based technique for translation, rotation and scaleinvariant image registration," IEEE Trans. Image Processing, 5(8), 1266–1271 (1996).
http://www.lira.dist.unige.it/teaching/SINA_10/slides-current/fourier-mellin-paper.pdf

[2] Manuel Guizar-Sicairos, Samuel T. Thurman, and James R. Fienup, "Efficient subpixel image registration algorithms," Opt. Lett. 33, 156-158 (2008).
https://pdfs.semanticscholar.org/b597/8b756bdcad061e3269eafaa69452a0c43e1b.pdf

[1]は平行移動、回転、スケーリングされた画像のピクセルを重ね合わせるアルゴリズム、[2]はサブピクセル精度で効率的に重ね合わせるアルゴリズムである。

Pythonのプログラムはアルゴリズム[2]のみ実装されたMatlabプログラムを参考にしている。
https://jp.mathworks.com/matlabcentral/fileexchange/18401-efficient-subpixel-image-registration-by-cross-correlation

3. ソースコードの解説

ソースコード↓
https://github.com/logicool-repo/fft-based_image_coregistration/blob/master/main.py

ソースコードは複数の関数とそれらを使用して実際に画像重ね合わせを行うmainから成る。

関数一覧
- fft_coreg_trans
- fft_coreg_LP
- dftregistration
- dftups
- FTpad
- logpolar_module
- main

3.1. main

まずは重ね合わせの流れを理解するためmainの中身を説明する。

slaveのシフト、回転、スケール量を設定

## set slave shift/rotate/scale values
trans_true = [2,-5]
angle_true = 2
scale_true = 1.002
mag_scale = 100

master, slave画像の読込と画像の縦横長、中心位置の取得

## load master and slave images
img_dir = 'lena/'
f = np.asarray(cv2.imread(img_dir+'lena512.png',0),dtype=np.float64)
g = np.asarray(cv2.imread(img_dir+'lena512.png',0),dtype=np.float64)
f = f[slice(512),slice(512)]
g = g[slice(512),slice(512)]

row = f.shape[0]; col = f.shape[1] # row & col size
hrow = int(row/2); hcol = int(col/2)
center = tuple(np.array(f.shape)/2)

slave画像の加工(平行移動、スケール、回転)

## translate slave
transMat = np.float32([[1,0,trans_true[0]],[0,1,trans_true[1]]])
g = cv2.warpAffine(g,transMat,(col,row))

## scale slave
g_tmp = cv2.resize(g,None,fx=scale_true,fy=scale_true, interpolation = cv2.INTER_LANCZOS4)
row_pad = int(g_tmp.shape[0]/2 - 512/2)
col_pad = int(g_tmp.shape[1]/2 - 512/2)

g = g*0
if scale_true < 1.0:
    row_tmp = g_tmp.shape[0]; col_tmp = g_tmp.shape[1] # row & col size
    hrow_tmp = np.floor(row_tmp/2); hcol_tmp = np.floor(col_tmp/2)
    if row_tmp % 2 == 0:
        row_slice = slice(int(center[0]-hrow_tmp),int(center[0]+hrow_tmp))
        col_slice = slice(int(center[1]-hcol_tmp),int(center[1]+hcol_tmp))
    else:
        row_slice = slice(int(center[0]-hrow_tmp),int(center[0]+hrow_tmp+1))
        col_slice = slice(int(center[1]-hcol_tmp),int(center[1]+hcol_tmp+1))
    g[row_slice,col_slice] = g_tmp
else:
    row_tmp = g_tmp.shape[0]; col_tmp = g_tmp.shape[1] # row & col size
    if row_tmp % 2 == 0:
        center_tmp = np.array(g_tmp.shape)/2
        row_slice = slice(int(center_tmp[0]-hrow),int(center_tmp[0]+hrow))
        col_slice = slice(int(center_tmp[1]-hcol),int(center_tmp[1]+hcol))
    else:
        center_tmp = np.floor(np.array(g_tmp.shape)/2)
        row_slice = slice(int(center_tmp[0]-hrow),int(center_tmp[0]+hrow))
        col_slice = slice(int(center_tmp[1]-hcol),int(center_tmp[1]+hcol))
    g = g_tmp[row_slice,col_slice]

## rotate slave
rotMat = cv2.getRotationMatrix2D(center, angle_true, 1.0)
g = cv2.warpAffine(g, rotMat, g.shape, flags=cv2.INTER_LANCZOS4)

master, slaveを周波数領域でLog-Polar変換

## Fourier log-magnitude spectra mapping in Log-Polar plane 
FLP, GLP = logpolar_module(f,g,mag_scale)

回転とスケール量を同時に推定

## estimate angle & scale
row_shift, col_shift, peak_map = fft_coreg_LP(FLP,GLP)
angle_est = - row_shift/(hrow) * 180
scale_est = 1.0 - col_shift/mag_scale

slaveを回転

## rotate slave
rotMat = cv2.getRotationMatrix2D(center, angle_est, 1.0)
g_coreg = cv2.warpAffine(g, rotMat, g.shape, flags=cv2.INTER_LANCZOS4)

続いてslaveをスケーリング、残るは平行移動の補正のみになる

## scale slave
g_coreg_tmp = cv2.resize(g_coreg,None,fx=scale_est,fy=scale_est, interpolation = cv2.INTER_LANCZOS4)
row_coreg_tmp = g_coreg_tmp.shape[0]; col_coreg_tmp = g_coreg_tmp.shape[1]
g_coreg = np.zeros((row,col))
if row_coreg_tmp == row:
    g_coreg = g_coreg_tmp
elif row_coreg_tmp > row:
    g_coreg = g_coreg_tmp[slice(row),slice(col)]
else:
    g_coreg[slice(row_coreg_tmp),slice(col_coreg_tmp)] = g_coreg_tmp

slaveの平行移動量の推定、slaveの平行移動

## estimate translation & translate slave
row_shift, col_shift, peak_map, g_coreg = fft_coreg_trans(f,g_coreg)

(後は画像プロットのコードなので省略)

3.2. fft_coreg_trans, fft_coreg_LP

fft_coreg_transは平行移動量を推定、fft_coreg_LPは回転、スケール量を推定する。これらの関数の中身はほぼ共通なのでfft_coreg_transだけ説明する。

2次元のFFTで画像を周波数領域に変換

## fft2
master_fd = fft2(master)
slave_fd = fft2(slave)

全画素の振幅を1にする(※画像のずれ=位相シフトだけが知りたいので)

## normalization
master_nfd = master_fd/np.abs(master_fd)
slave_nfd = slave_fd/np.abs(slave_fd)

サブピクセル単位のズレ量を推定

## shift estimation
usfac = 100
output, Nc, Nr, peak_map = dftregistration(master_nfd,slave_nfd,usfac)

nr, nc = slave.shape
diffphase = output[1]
row_shift = output[2]
col_shift = output[3]

推定量に基づいてslaveを平行移動する(※fft_coreg_LPの場合はmainで回転とスケールする)

## coregistration
slave_fd_crg = slave_fd*np.exp(1j*2*np.pi*(-row_shift*Nr/nr-col_shift*Nc/nc))*np.exp(1j*diffphase)
slave_crg = ifft2(slave_fd_crg)
slave_crg = np.abs(slave_crg)

3.2. dftregistration

これが重ね合わせアルゴリズムのアルゴリズムの肝となる関数である。

3.3. dftups

3.4. FTpad

3.5. logpolar_module

回転とスケール量を推定するために周波数領域のmasterとslave画像のLog-Polar変換結果を出力するだけの関数。slaveの平行移動、回転、スケールに伴い、画像の外縁部分は重ね合わせられないので窓関数を適用してから2次元FFTを適用する。続いて、周波数領域のデータに対してhighpass filterを適用する。最後に、Log-Polar変換を施す。ここではOpenCVの関数を使っているが、特に特殊な処理ではないのでライブラリを使ってしまって全く問題ない。

row = f.shape[0]; col = f.shape[1] # row & col size
hrow = int(row/2); hcol = int(col/2)

## hanning window
hy = np.hanning(row)
hx = np.hanning(col)
hw = hy.reshape(row, 1) * hx.reshape(1, col)
f = f * hw
g = g * hw

# fft
F = fftshift(fft2(f))
G = fftshift(fft2(g))

# highpass filter
X1 = np.cos(np.pi*(np.arange(row)/row-0.5))
X2 = np.cos(np.pi*(np.arange(col)/col-0.5))
X1 = np.reshape(X1,(row,1))
X2 = np.reshape(X2,(1,col))
X1 = np.tile(X1,(1,col))
X2 = np.tile(X2,(row,1))
X = X1*X2
H = (1.0-X)*(2.0-X)
F = H * F
G = H * G

## Log-Polar transform
F = np.abs(F)
G = np.abs(G)
FLP = cv2.logPolar(F, (F.shape[0]/2, F.shape[1]/2), mag_scale, cv2.INTER_LANCZOS4)
GLP = cv2.logPolar(G, (G.shape[0]/2, G.shape[1]/2), mag_scale, cv2.INTER_LANCZOS4)

## roll and slice
FLP = np.roll(FLP,int(hcol),axis=1)
GLP = np.roll(GLP,int(hcol),axis=1)
FLP = FLP[slice(int(hrow)),slice(int(hcol))]
GLP = GLP[slice(int(hrow)),slice(int(hcol))]

変換結果FLP、GLPの右半分だけを出力しているのは、Log-Polar画像のシフト量計算で全領域を使ってしまうと支障があるためである。

Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away