2-D cross-correlation - MATLAB xcorr2 (2024)

Create two matrices, M1 and 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 is 5-by-5 and M2 is 3-by-3, so their cross-correlation has size (5+3-1)-by-(5+3-1), or 7-by-7. In terms of lags, the resulting matrix is

$$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).$$

As an example, compute the element $$c_{0,2}$$ (or C(3,5) in MATLAB®, since M2 is 3-by-3). Line up the two matrices so their (1,1) elements coincide. This placement corresponds to $$c_{0,0}$$. To find $$c_{0,2}$$, slide M2 two rows to the right.

Now M2 is on top of the matrix M1(1:3,3:5). Compute the element-by-element products and sum them. The answer should be

$$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

Verify the result using xcorr2.

D = xcorr2(M1,M2);DD = D(0+r2,2+c2)

DD = 585

Given a matrix $$\mathcal{X}$$ of size $$M\times N$$ and a matrix $$\mathcal{H}$$ of size $$P\times Q$$, their two-dimensional cross-correlation, $$\mathcal{C}=\mathcal{X}\star \mathcal{H}$$, is a matrix of size $$(M+P-1)\times (N+Q-1)$$ with elements

$$\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$$ is the trace and the dagger denotes Hermitian conjugation. The matrices $$\underset{}{\overset{\sim }{\mathcal{X}}}$$ and $${\underset{}{\overset{\sim }{\mathcal{H}}}}_{kl}$$ have size $$(M+2(P-1))\times (N+2(Q-1))$$ and nonzero elements given by

$$\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}$$

and

$${\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}$$

Calling xcorr2 is equivalent to this procedure for general complex matrices of arbitrary size.

Create two complex matrices, $$\mathcal{X}$$ of size $$7\times 22$$ and $$\mathcal{H}$$ of size $$6\times 17$$.

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;

Initialize $$\underset{}{\overset{\sim }{\mathcal{X}}}$$ and $$\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]);

Compute the elements of $$\mathcal{C}$$ by looping over $$k$$ and $$l$$. Reset $${\underset{}{\overset{\sim }{\mathcal{H}}}}_{kl}$$ to zero at each step. Save time and memory by summing element products instead of multiplying and taking the trace.

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))); endendmax(max(abs(C-xcorr2(X,H))))

ans = 1.5139e-14

The answer coincides to machine precision with the output of xcorr2.

Use cross-correlation to find where a section of an image fits in the whole. Cross-correlation enables you to find the regions in which two signals most resemble each other. For two-dimensional signals, like images, use xcorr2.

Load a black-and-white test image into the workspace. Display it with imagesc.

load durerimg = X;White = max(max(img));imagesc(img)axis image offcolormap graytitle('Original')

Select a rectangular section of the image. Display the larger image with the section missing.

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 offcolormap graytitle('Image')subplot(1,2,2)imagesc(kumg)axis image offcolormap graytitle('Section')

Use xcorr2 to find where the small image fits in the larger image. Subtract the mean value so that there are roughly equal numbers of negative and positive values.

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

The maximum of the cross-correlation corresponds to the estimated location of the lower-right corner of the section. Use ind2sub to convert the one-dimensional location of the maximum to two-dimensional coordinates.

[ssr,snd] = max(crr(:));[ij,ji] = ind2sub(size(crr),snd);figureplot(crr(:))title('Cross-Correlation')hold onplot(snd,ssr,'or')hold offtext(snd*1.05,ssr,'Maximum')

Place the smaller image inside the larger image. Rotate the smaller image to comply with the convention that MATLAB® uses to display images. Draw a rectangle around it.

img(ij:-1:ij-size(Sect,1)+1,ji:-1:ji-size(Sect,2)+1) = rot90(Sect,2);imagesc(img)axis image offcolormap graytitle('Reconstructed')hold onplot([y y Y Y y],[x X X x x],'r')hold off

Shift a template by a known amount and recover the shift using cross-correlation.

Create a template in an 11-by-11 matrix. Create a 22-by-22 matrix and shift the original template by 8 along the row dimension and 6 along the column dimension.

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;

Plot the original and shifted templates.

imagesc(offsetTemplate)colormap grayhold onimagesc(template)axis equal

Cross-correlate the two matrices and find the maximum absolute value of the cross-correlation. Use the position of the maximum absolute value to determine the shift in the template. Check the result against the known shift.

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

The shift obtained from the cross-correlation equals the known template shift in the row and column dimensions.

This example requires Parallel Computing Toolbox™ software. Refer to GPU Computing Requirements (Parallel Computing Toolbox) to see what GPUs are supported.

Shift a template by a known amount and recover the shift using cross-correlation.

Create a template in an 11-by-11 matrix. Create a 22-by-22 matrix and shift the original template by 8 along the row dimension and 6 along the column dimension.

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;

Put the original and shifted template matrices on your GPU using gpuArray objects.

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

Compute the cross-correlation on the GPU.

cc = xcorr2(offsetTemplate,template);

Return the result to the MATLAB® workspace using gather. Use the maximum absolute value of the cross-correlation to determine the shift, and compare the result with the known shift.

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

c — 2-D cross-correlation or autocorrelation matrix
matrix

2-D cross-correlation or autocorrelation matrix, returned as a matrix.