xcorr2

2-D cross-correlation

Syntax

C = xcorr2(A,B)
A = xcorr2(A)
C = xcorr2(gpuArrayA,gpuArrayB)

Description

C = xcorr2(A,B) returns the cross-correlation of matrices A and B with no scaling. xcorr2 is the two-dimensional version of xcorr.

A = xcorr2(A) is the autocorrelation matrix of input matrix A. It is identical to xcorr2(A,A).

C = xcorr2(gpuArrayA,gpuArrayB) returns the cross-correlation of two matrices of class gpuArray. The output cross-correlation matrix, C, is also a gpuArray object. See Establish Arrays on a GPU for details on gpuArray objects. Using xcorr2 with gpuArray objects requires Parallel Computing Toolbox™ software and a CUDA-enabled NVIDIA GPU with compute capability 1.3 or above. See http://www.mathworks.com/products/parallel-computing/requirements.html for details. See GPU Acceleration for Cross-Correlation Matrix Computation for an example of using a GPU to compute the cross-correlation.

2-D Cross-Correlation

The 2-D cross-correlation of an M-by-N matrix X and a P-by-Q matrix H is a matrix C of size M+P–1 by N+Q–1 given by

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

where the bar over H denotes complex conjugation.

The output matrix, C(k,l), has negative and positive row and column indices. A negative row index corresponds to an upward shift of the rows of H. A negative column index corresponds to a leftward shift of the columns of H. A positive row index corresponds to a downward shift of the rows of H. A positive column index corresponds to a rightward shift of the columns. To cast the indices in MATLAB® form, simply add the size of H: the element C(k,l) corresponds to C(k+P,l+Q) in the workspace.

For example, consider the following 2-D cross-correlation:

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

The C(1,1) element in the output corresponds to C(1–3,1–2) = C(–2,–1) in the defining equation, which uses zero-based indexing. The C(1,1) element is computed by shifting H two rows upward and one column to the left. Accordingly, the only product in the cross-correlation sum is X(1,1)*H(3,2) = 6. Using the defining equation, you obtain

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

with all other terms in the double sum equal to zero.

Examples

Output Matrix Size

If matrix I1 has dimensions (4,3) and matrix I2 has dimensions (2,2), the following equations determine the number of rows and columns of the output matrix:

Cfullrows=I1rows+I2rows1=4+21=5Cfullcolumns=I1columns+I2columns1=3+21=4

The resulting matrix is

Cfull=[c00c01c02c03c10c11c12c13c20c21c22c23c30c31c32c33c40c41c42c43]

Computing a Specific Element

Cvalidcolumns=I1columnsI2columns+1=2

In cross-correlation, the value of an output element is computed as a weighted sum of neighboring elements. For example, suppose the first input matrix represents an image and is defined as

I1 = [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]

The second input matrix also represents an image and is defined as

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

The following figure shows how to compute the (2,4) output element (zero-based indexing) using these steps:

  1. Slide the center element of I2 so that lies on top of the (1,3) element of I1.

  2. Multiply each weight in I2 by the element of I1 underneath.

  3. Sum the individual products from step 2.

The (2,4) output element from the cross-correlation is

18+81+156+73+145+167+134+209+222=585

The normalized cross-correlation of the (2,4) output element is

585/sqrt(sum(dot(I1p,I1p))*sum(dot(I2,I2))) = 0.8070

where I1p = [1 8 15; 7 14 16; 13 20 22].

Recovery of Template Shift with Cross-Correlation

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 = .2*ones(11);
template(6,3:9) = .6;
template(3:9,6) = .6;
offsetTemplate = .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 gray;
hold on;
imagesc(template);

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)

The returned 1 indicates that the shift obtained the cross-correlation equals the known the template shift in both the row and column dimension.

GPU Acceleration for Cross-Correlation Matrix Computation

The following example requires Parallel Computing Toolbox software and a CUDA-enabled NVIDIA GPU with compute capability 1.3 or above. See http://www.mathworks.com/products/parallel-computing/requirements.html for details.

Repeat the example Recovery of Template Shift with Cross-Correlation. For convenience, the code to create the original and shifted templates is repeated.

template = .2*ones(11); 
template(6,3:9) = .6;   
template(3:9,6) = .6;
offsetTemplate = .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) 

See Also

| |

Was this topic helpful?