Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

xcorr2

説明

c = xcorr2(a,b) は、スケーリングなしで、行列 ab の相互相関を返します。xcorr2xcorr の 2 次元バージョンです。

c = xcorr2(a) は、入力行列 a の自己相関行列です。この構文は、xcorr2(a,a) と等価です。

すべて折りたたむ

2 つの行列 M1 および M2 を作成します。

M1 = [17 24  1  8 15;
      23  5  7 14 16;
       4  6 13 20 22;
      10 12 19 21  3;
      11 18 25  2  9];

M2 = [8 1 6;
      3 5 7;
      4 9 2];

M1 は 5 行 5 列で M2 は 3 行 3 列であることから、これらの相互相関のサイズは (5+3-1) 行 (5+3-1) 列、つまり 7 行 7 列となります。ラグについては、結果の行列は以下のようになります。

C=(c-2,-2c-2,-1c-2,0c-2,1c-2,2c-2,3c-2,4c-1,-2c-1,-1c-1,0c-1,1c-1,2c-1,3c-1,4c0,-2c0,-1c0,0c0,1c0,2c0,3c0,4c1,-2c1,-1c1,0c1,1c1,2c1,3c1,4c2,-2c2,-1c2,0c2,1c2,2c2,3c2,4c3,-2c3,-1c3,0c3,1c3,2c3,3c3,4c4,-2c4,-1c4,0c4,1c4,2c4,3c4,4).

例として、要素 c0,2 (または M2 が 3 行 3 列のため、MATLAB® における C(3,5)) を計算します。(1,1) の要素が一致するように 2 つの行列を整列させます。この配置は c0,0 に対応します。c0,2 を見つけるには、M2 の 2 つの行を右にスライドさせます。

これで M2 が行列 M1(1:3,3:5) の上に重なりました。要素ごとの積を計算してそれらを加算します。答えは以下のようになります。

1×8+7×3+13×4+8×1+14×5+20×9+15×6+16×7+22×2=585.

[r2,c2] = size(M2);

CC = sum(sum(M1(0+(1:r2),2+(1:c2)).*M2))
CC = 585

xcorr2 を使用して結果を確認します。

D = xcorr2(M1,M2);

DD = D(0+r2,2+c2)
DD = 585

サイズ M×N の行列 X およびサイズ P×Q の行列 H が与えられた場合、それらの 2 次元相互相関 C=XH は以下の要素を持つサイズ (M+P-1)×(N+Q-1) の行列になります。

C(k,l)=Tr{XHkl}1kM+P-1,1lN+Q-1.

Tr はトレースで、ダガーはエルミート共役を意味します。行列 X および Hkl のサイズは (M+2(P-1))×(N+2(Q-1)) であり、その非ゼロ要素は以下で与えられます。

X(m,n)=X(m-P+1,n-Q+1),PmM+P-1,QnN+Q-1

Hkl(p,q)=H(p-k+1,q-l+1),kpP+k-1,lqQ+l-1.

xcorr2 の呼び出しは任意のサイズの通常の複素行列についてのこの処理手順と等価です。

サイズ 7×22X およびサイズ 6×17H という 2 つの複素行列を作成します。

X = randn([7 22])+1j*randn([7 22]);
H = randn([6 17])+1j*randn([6 17]);

[M,N] = size(X);
m = 1:M;
n = 1:N;

[P,Q] = size(H);
p = 1:P;
q = 1:Q;

XC を初期化します。

Xt = zeros([M+2*(P-1) N+2*(Q-1)]);
Xt(m+P-1,n+Q-1) = X;
C = zeros([M+P-1 N+Q-1]);

k および l をループ処理することによって C の要素を計算します。各ステップで Hkl をゼロにリセットします。乗算してトレースを取得する代わりに要素の積を合計することによって、時間とメモリを節約します。

for k = 1:M+P-1
    for l = 1:N+Q-1
        Hkl = zeros([M+2*(P-1) N+2*(Q-1)]);
        Hkl(p+k-1,q+l-1) = H;
        C(k,l) = sum(sum(Xt.*conj(Hkl)));
    end
end

max(max(abs(C-xcorr2(X,H))))
ans = 1.5139e-14

答えは、xcorr2 の出力のマシンの精度と一致しています。

相互相関を使用して、イメージの一部が全体のどこに収まるかを見つけます。相互相関により、含まれる 2 つの信号が互いに最もよく似ている領域を検出できます。イメージのような 2 次元信号に対しては xcorr2 を使用します。

白黒のテスト イメージをワークスペースに読み込みます。imagesc でこれを表示します。

load durer
img = X;
White = max(max(img));

imagesc(img)
axis image off
colormap gray
title('Original')

Figure contains an axes object. The axes object with title Original contains an object of type image.

イメージの一部分を四角形に選択します。この部分を取り除いた大きい方のイメージを表示します。

x = 435;
X = 535;
szx = x:X;

y = 62;
Y = 182;
szy = y:Y;

Sect = img(szx,szy);

kimg = img;
kimg(szx,szy) = White;

kumg = White*ones(size(img));
kumg(szx,szy) = Sect;

subplot(1,2,1)
imagesc(kimg)
axis image off
colormap gray
title('Image')

subplot(1,2,2)
imagesc(kumg)
axis image off
colormap gray
title('Section')

Figure contains 2 axes objects. Axes object 1 with title Image contains an object of type image. Axes object 2 with title Section contains an object of type image.

xcorr2 を使用して、小さい方のイメージが大きい方のイメージのどこに当てはまるかを見つけます。正と負の値の数がほぼ等しくなるように平均値を差し引きます。

nimg = img-mean(mean(img));
nSec = nimg(szx,szy);

crr = xcorr2(nimg,nSec);

相互相関の最大値は小さい方のイメージの右下隅の推定位置に相当します。ind2sub を使用して、最大値の 1 次元的位置を 2 次元座標に変換します。

[ssr,snd] = max(crr(:));
[ij,ji] = ind2sub(size(crr),snd);

figure
plot(crr(:))
title('Cross-Correlation')
hold on
plot(snd,ssr,'or')
hold off
text(snd*1.05,ssr,'Maximum')

Figure contains an axes object. The axes object with title Cross-Correlation contains 3 objects of type line, text. One or more of the lines displays its values using only markers

大きい方のイメージ内に小さい方のイメージを配置します。MATLAB® がイメージ表示で使用する方法に従い、小さい方のイメージを回転させます。小さい方のイメージの周囲に四角形を描きます。

img(ij:-1:ij-size(Sect,1)+1,ji:-1:ji-size(Sect,2)+1) = rot90(Sect,2);

imagesc(img)
axis image off
colormap gray
title('Reconstructed')
hold on
plot([y y Y Y y],[x X X x x],'r')
hold off

Figure contains an axes object. The axes object with title Reconstructed contains 2 objects of type image, line.

既知の量だけテンプレートをシフトし、相互相関を使用してシフトを回復します。

11 行 11 列の行列をもつテンプレートを作成します。22 行 22 列の行列を作成してその行の次元に沿って 8、列の次元に沿って 6 だけ元のテンプレートをシフトします。

template = 0.2*ones(11);
template(6,3:9) = 0.6;
template(3:9,6) = 0.6;
offsetTemplate = 0.2*ones(22);
offset = [8 6];
offsetTemplate((1:size(template,1))+offset(1), ...
    (1:size(template,2))+offset(2)) = template;

元の行列とシフトした行列をプロットします。

imagesc(offsetTemplate)
colormap gray
hold on
imagesc(template)
axis equal

Figure contains an axes object. The axes object contains 2 objects of type image.

2 つの行列を相互相関させ、相互相関の最大絶対値を見つけます。最大絶対値の位置を使用して、テンプレート内のシフトを決定します。既知のシフトに対して結果を確認します。

cc = xcorr2(offsetTemplate,template);
[max_cc, imax] = max(abs(cc(:)));
[ypeak, xpeak] = ind2sub(size(cc),imax(1));
corr_offset = [(ypeak-size(template,1)) (xpeak-size(template,2))];

isequal(corr_offset,offset)
ans = logical
   1

相互相関から得られたシフトは、行と列の次元において既知のテンプレートのシフトと等価です。

この例では、Parallel Computing Toolbox™ ソフトウェアが必要です。どの GPU がサポートされているかについては、GPU 計算の要件 (Parallel Computing Toolbox)を参照してください。

既知の量だけテンプレートをシフトし、相互相関を使用してシフトを回復します。

11 行 11 列の行列をもつテンプレートを作成します。22 行 22 列の行列を作成してその行の次元に沿って 8、列の次元に沿って 6 だけ元のテンプレートをシフトします。

template = 0.2*ones(11); 
template(6,3:9) = 0.6;   
template(3:9,6) = 0.6;
offsetTemplate = 0.2*ones(22); 
offset = [8 6];
offsetTemplate((1:size(template,1))+offset(1), ...
    (1:size(template,2))+offset(2)) = template;

gpuArray オブジェクトを使用して、GPU 上に元のテンプレートとシフトされたテンプレートの行列を配置します。

template = gpuArray(template);
offsetTemplate = gpuArray(offsetTemplate);

GPU 上での相互相関を計算します。

cc = xcorr2(offsetTemplate,template);

gather を使用して、結果を MATLAB® ワークスペースに返します。相互相関の最大絶対値を使用してシフトを決定し、結果を既知のシフトと比較します。

cc = gather(cc);
[max_cc,imax] = max(abs(cc(:)));
[ypeak,xpeak] = ind2sub(size(cc),imax(1));
corr_offset = [(ypeak-size(template,1)) (xpeak-size(template,2))];
isequal(corr_offset,offset)
ans = logical
   1

入力引数

すべて折りたたむ

入力配列。行列として指定します。

例: sin(2*pi*(0:9)'/10)*sin(2*pi*(0:13)/20) は、2 次元の正弦波表面を指定します。

データ型: single | double
複素数のサポート: あり

出力引数

すべて折りたたむ

2 次元の相互相関行列または自己相関行列。行列として返されます。

詳細

すべて折りたたむ

2 次元相互相関

M 行 N 列の行列 X と P 行 Q 列の行列 H の 2 次元相互相関は、サイズが M+P–1 行 N+Q–1 列の行列 C となります。この要素は次で求められます。

C(k,l)=m=0M1n=0N1X(m,n)H¯(mk,nl),(P1)kM1,(Q1)lN1,

H の上の横線は複素共役を表します。

出力される行列 C(k,l) には、負と正の行インデックスと列インデックスがあります。

  • 負の行インデックスは、H の行の上方向のシフトに対応します。

  • 負の列インデックスは、H の列の左方向のシフトに対応します。

  • 正の行インデックスは、H の行の下方向のシフトに対応します。

  • 正の列インデックスは、H の列の右方向のシフトに対応します。

インデックスを MATLAB® 形式にキャストするには、H のサイズを加算します。要素 C(k,l) は、ワークスペースの C(k+P,l+Q) に対応しています。

たとえば、以下の 2 次元相互相関について考えます。

X = ones(2,3);
H = [1 2; 3 4; 5 6];  % H is 3 by 2
C = xcorr2(X,H)
C =
     6    11    11     5
    10    18    18     8
     6    10    10     4
     2     3     3     1

出力における要素 C(1,1) は定義式の C(1–3,1–2) = C(–2,–1) に対応し、ゼロベースのインデックスを使用します。C(1,1) 要素を計算するには、H を 2 行上、1 列左にシフトします。その結果、相互相関の和における積は唯一 X(1,1)*H(3,2) = 6 となります。定義式を使用すると、

C(2,1)=m=01n=02X(m,n)H¯(m+2,n+1)=X(0,0)H¯(2,1)=1×6=6,

が得られ、二重和の他のすべての項はゼロになります。

拡張機能

C/C++ コード生成
MATLAB® Coder™ を使用して C および C++ コードを生成します。

バージョン履歴

R2006a より前に導入

参考

| |