画像の平行移動位置あわせ（差の二乗を使用）の例

前準備

Octave のインストールが済んでいること．

必見 Web ページ: http://www.csse.uwa.edu.au/~pk/Research/MatlabFns/

必見 Web ページ: http://www.eecs.berkeley.edu/Research/Projects/CS/vision/bsds/

* Octave のプロンプトを変えたいときは，次のように操作する．

PS1('> ')

* Octave のインストールによっては，Octave の起動時に毎回次の操作を行う必要があるかもしれない

pkg load image

プログラムのソースコード

# 2つの濃淡画像について，差の二乗の総和が最小になるように平行移動位置合わせ
#   画像を平行移動し，2つの画像の重なる部分の全画素についての，差の二乗の総和を求め，それを画素数で割る．それが最小になるような平行移動量を求めたい．

function err = compare_mat( A, B, dy, dx )
  # 引数の並びは dy が先で dx が後
  # dx > 0 かつ dy > 0 のときは，
  # A の右下部分と B の左上部分との比較.
  # 「A の右下部分」というのは，A の上 dy 行，左 dy 列を取り除いた残り
  # 「B の左下部分」と「A の右下部分」が同じサイズになるように
  #  A, B から部分を取り出して比較に使う．
  # dx, dy が負のときは「A の右下部分」が左上、左下、左上に変わる
  # 縦と横を知る
  heightA = size(A)(1);
  widthA = size(A)(2);
  heightB = size(B)(1);
  widthB = size(B)(2);

  if ( ( dx >= 0 ) & (dy >= 0) )
    # A 左上 から dx, dy 分は使わない．残りを AA へ.
    # dx = 0, dy = 0 ならば A と AA は等しい
    AA = A(1 + dy:heightA,          1 + dx:widthA); 
    # BB は 左上 から AA と同じ分だけ使う
    BB = B(1:size(AA)(1),           1:size(AA)(2));
  elseif ( ( dx >= 0 ) & (dy < 0) )
    AA = A(1:heightA + dy,          1 + dx:widthA); 
    BB = B(1 - dy:size(AA)(1) - dy, 1:size(AA)(2));
  elseif ( ( dx < 0 ) & (dy >= 0) )
    AA = A(1 + dy:heightA,          1:widthA + dx); 
    BB = B(1:size(AA)(1),           1 - dx:size(AA)(2) - dx);
  elseif ( ( dx < 0 ) & (dy < 0) )
    AA = A(1:heightA + dy,          1:widthA + dx); 
    BB = B(1 - dy:size(AA)(1) - dy, 1 - dx:size(AA)(2) - dx);
  endif

# 各要素の差を取り，その2乗の総和を求め、AA のサイズで割る
  Z = AA - BB;
  err = sum( sum( Z .* Z ) ) / ( size(AA)(1) * size(AA)(2) );
  clear AA;
  clear BB;
  clear Z;
endfunction

function [x, y, val] = registration(gray1, gray2, y1, y2, x1, x2)
# 画像の平行移動位置あわせ（総当たりによる方法）
# x, y が移動量．val はそのときの「相違度」

# 2つの濃淡画像gray1, gray2 で，「相違度」が小さくなるように，画像をずらしたい．
# 関数 compare_mat を補助的に使い，
# y1, y2, x1, x2 は移動の範囲．
# この範囲の中をすべて試し，最も「相違度」が小さいものを評価結果 [x, y, val] とする．

# 結果の格納用の配列を作る
  d = ones(y2 - y1 + 1, x2 - x1 + 1) * 256;

  # rangeX, rangeY はずらす量
  # 単純な総当たり
  rangeX = [x1:x2];
  rangeY = [y1:y2];
  for i = [1:size(rangeX)(2)]
    for j = [1:size(rangeY)(2)]
      # 行列 d には相違度の値を格納
      d(j, i) = compare_mat( gray1, gray2, rangeY(j), rangeX(i) );
    endfor
  endfor
  # 最も小さい値を持つのはどこか？
  [val_vec, index_vec] = min(d);
  [val, ix] = min(val_vec);
  iy = index_vec(ix);
  x = rangeX(ix);
  y = rangeY(iy);
  clear d;
  clear rangeX;
  clear rangeY;
endfunction

function [x, y, val] = registration_random(gray1, gray2, Nsamples, y1, y2, x1, x2)
# 画像の平行移動位置あわせ（ランダムサンプリングによる方法）
# x, y が移動量．val はそのときの「相違度」

# 2つの濃淡画像gray1, gray2 で，「相違度」が小さくなるように，画像をずらしたい．
# 関数 compare_mat を補助的に使い，
# ランダムに試してみる．
# Nsample はサンプル数，y1, y2, x1, x2 は移動の範囲．
# このサンプル数の試行の中で最も「相違度」が小さいものを評価結果 [x, y, val] とする．
# 一様乱数 Nsamples 個ずつ．rx, ry は列ベクトル(N行1列のベクトル)
  rx = ceil( rand(Nsamples,1) * ((x2 - x1) + 1) ) + (x1 - 1);
  ry = ceil( rand(Nsamples,1) * ((y2 - y1) + 1) ) + (y1 - 1);

  # 試す
  for i = [1:Nsamples]
    similar(i) = compare_mat( gray1, gray2, ry(i), rx(i) );
  endfor

  [val, index] = min(similar);
  x = rx(index(1));
  y = ry(index(1));
endfunction

function [x, y, v] = registration_pyramid(gray1, gray2, N )
# find the minimum using pyramid way

# まず「小さい」画像を作って，それで位置あわせを行うことにしたい
#   その小さい画像のサイズを求める．N より小さくならないようにする
MI = min( [size(gray1)(1) size(gray1)(2) size(gray2)(1) size(gray2)(2)] );
#   REPEAT は，「1/2 に縮小する」処理を何回行うと「小さい」画像になるか
REPEAT = floor( log2(MI) - log2(N) )
#   「小さい」画像 A, B を作る
A = imresize( gray1, 1/(2^REPEAT), 'cubic' );
B = imresize( gray2, 1/(2^REPEAT), 'cubic' );

# せめて A, B の上下 1/3 が重なり，左右 1/3 が重なるように
dy = floor( min( [size(A)(1) * 0.67 size(B)(1) * 0.67] ) );
dx = floor( min( [size(A)(2) * 0.67 size(B)(2) * 0.67] ) );
[rx, ry, val] = registration(A, B, -dy, dy, -dx, dx);
clear A;
clear B;

for i = (REPEAT-1):-1:0
  disp "."
  # 前回よりも縦横２倍大きい画像
  A = imresize( gray1, 1/(2^i), 'cubic' );
  B = imresize( gray2, 1/(2^i), 'cubic' );
  # それほど広い範囲を探す必要は無い．２画素分でいいだろう
  [new_rx, new_ry, val] = registration(A, B, ry * 2 - 2, ry * 2 + 2, rx * 2 - 2, ry * 2 + 2); 
  rx = new_rx;
  ry = new_ry;
  clear A;
  clear B;
endfor

x = rx;
y = ry;
v = val;
endfunction

# 画像ファイル 2つを読み込み
rgb1 = imread("/tmp/441.png");
rgb2 = imread("/tmp/442.png");
gray1 = rgb2gray( rgb1 );
gray2 = rgb2gray( rgb2 );

[x y v] = registration_pyramid(gray1, gray2, 64);
x
y
v


# For debug and test of functions
# compare_mat( gray1, gray2, 0, 1 );
# # [rx, ry, val]    = registration(gray1, gray2, 0, 50, 0, 50);
# rx
# ry
# val
# # [rx2, ry2, val2] = registration_random(gray1, gray2, 1000, 0, 50, 0, 50);
# rx
# ry
# val2

実行方法と実行結果の例

次の 2つの画像の平行移動位置あわせを行っている例です．（この画像は lena 標準画像から切り抜いて作ったもの）． 1つ目の画像の変異量（左が正で、右が負で、上が正で、下が負）を求めます．

この種のソフトは，探せば同様のが見つかるとも思いますが， Octave の中に組み込んで使う（つまり，平行移動位置あわせの結果をさらにいろいろな用途に使う）とか，平行移動位置あわせのやり方を細かく調整することをしたくて，Octave のソースコードプログラムを自前で書くことにした．

* この 2つの画像のファイルは，前準備として /tmp に保存しておくこと．

《１つ目》

《２つ目》

（オプション）Cygwin からリモートログインする場合
startx xhost + ssh -X username@ipaddress
ソースコードをファイルとして保存
Octave を起動

Octave の source コマンドを使って実行

cd ＜ソースコードのファイルを置いたディレクトリ＞
source "＜ソースコードのファイル名＞"

実行が終わるまで待つ