xcorr2

2 次元相互相関

構文

c = xcorr2(a,b)

c = xcorr2(a)

説明

c = xcorr2(a,b) は、スケーリングなしで、行列 a と b の相互相関を返します。xcorr2 は xcorr の 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=\left(\begin{array}{ccccccc}{c}_{-2,-2}& c_{-2,-1}& c_{-2,0}& c_{-2,1}& c_{-2,2}& c_{-2,3}& c_{-2,4}\\ c_{-1,-2}& c_{-1,-1}& c_{-1,0}& c_{-1,1}& c_{-1,2}& c_{-1,3}& c_{-1,4}\\ c_{0,-2}& c_{0,-1}& c_{0,0}& c_{0,1}& c_{0,2}& c_{0,3}& c_{0,4}\\ c_{1,-2}& c_{1,-1}& c_{1,0}& c_{1,1}& c_{1,2}& c_{1,3}& c_{1,4}\\ c_{2,-2}& c_{2,-1}& c_{2,0}& c_{2,1}& c_{2,2}& c_{2,3}& c_{2,4}\\ c_{3,-2}& c_{3,-1}& c_{3,0}& c_{3,1}& c_{3,2}& c_{3,3}& c_{3,4}\\ c_{4,-2}& c_{4,-1}& c_{4,0}& c_{4,1}& c_{4,2}& c_{4,3}& c_{4,4}\end{array}\right).$$

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

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

$$1\times 8+7\times 3+13\times 4+8\times 1+14\times 5+20\times 9+15\times 6+16\times 7+22\times 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

任意の複素行列の 2 次元相互相関

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

$$\mathcal{C}(k,l)=Tr\left\{\underset{}{\overset{\sim }{\mathcal{X}}}{\underset{}{\overset{\sim }{\mathcal{H}}}}_{kl}^{†}\right\}\begin{array}{c}1\le k\le M+P-1,\\ 1\le l\le N+Q-1.\end{array}$$

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

$$\underset{}{\overset{～}{\mathcal{X}}}(m,n)=\mathcal{X}(m-P+1,n-Q+1),\begin{array}{c}P\le m\le M+P-1,\\ Q\le n\le N+Q-1\end{array}$$

と

$${\underset{}{\overset{～}{\mathcal{H}}}}_{kl}(p,q)=\mathcal{H}(p-k+1,q-l+1),\begin{array}{c}k\le p\le P+k-1,\\ l\le q\le Q+l-1.\end{array}$$

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

サイズ $$7\times 22$$ の $$\mathcal{X}$$ およびサイズ $$6\times 17$$ の $$\mathcal{H}$$ という 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;

$$\underset{}{\overset{\sim }{\mathcal{X}}}$$ と $$\mathcal{C}$$ を初期化します。

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$$ をループ処理することによって $$\mathcal{C}$$ の要素を計算します。各ステップで $${\underset{}{\overset{\sim }{\mathcal{H}}}}_{kl}$$ をゼロにリセットします。乗算してトレースを取得する代わりに要素の積を合計することによって、時間とメモリを節約します。

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 つの信号が互いに最もよく似ている領域を検出できます。イメージのような 2 次元信号に対しては xcorr2 を使用します。

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

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

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

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

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')

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')

大きい方のイメージ内に小さい方のイメージを配置します。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

相互相関を使用したテンプレート シフトの回復

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

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

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

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

相互相関行列計算のための GPU 高速化

この例では、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

入力引数

a, b — 入力配列

行列

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

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

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

出力引数

c — 2 次元相互相関行列または自己相関行列

行列

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

詳細

2 次元相互相関

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

$$C(k,l)={\displaystyle \sum _{m=0}^{M-1}{\displaystyle \sum _{n=0}^{N-1}X(m,n)\overline{H}(m-k,n-l),}}\begin{array}{c}-(P-1)\le k\le M-1,\\ -(Q-1)\le l\le N-1,\end{array}$$

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)={\displaystyle \sum _{m=0}^1{\displaystyle \sum _{n=0}^{2}X(m,n)\overline{H}(m+2,n+1)}}=X(0,0)\overline{H}(2,1)=1\times 6=6,$$

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

拡張機能

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

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

この関数は、GPU 配列を完全にサポートします。詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2006a より前に導入

参考

conv2 | filter2 | xcorr