トップページ -> 実践知識 -> 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 "."
  # 前回よりも縦横2倍大きい画像
  A = imresize( gray1, 1/(2^i), 'cubic' );
  B = imresize( gray2, 1/(2^i), 'cubic' );
  # それほど広い範囲を探す必要は無い.2画素分でいいだろう
  [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 に保存しておくこと.

《1つ目》

《2つ目》

  1. (オプション)Cygwin からリモートログインする場合
    startx
    xhost +
    ssh -X username@ipaddress
    

  2. ソースコードをファイルとして保存

  3. Octave を起動

  4. Octave の source コマンドを使って実行
    cd <ソースコードのファイルを置いたディレクトリ>
    source "<ソースコードのファイル名>"
    

  5. 実行が終わるまで待つ