SignalProcessing - Maple Programming Help

Home : Support : Online Help : Science and Engineering : Signal Processing : Convolution and Correlation Computations : SignalProcessing/CrossCorrelation2D

SignalProcessing

 CrossCorrelation2D
 calculate the cross-correlation between a pair of one- or two-dimensional rtables

 Calling Sequence CrossCorrelation2D( A, B, options )

Parameters

 A, B - two one- or two-dimensional rtables container - (optional) container to compute and store the cross-correlation

Description

 • The CrossCorrelation2D command takes two 1-D or 2-D rtables and computes their 2-D cross-correlation.
 • The inputs are converted to Matrices of complex[8] datatype, and an error will be thrown if this is not possible. For this reason, it is most efficient for the inputs to already be Matrices having datatype complex[8].
 • The cross-correlation is computed via Discrete or Fast Fourier Transforms (DFTs or FFTs), but the result is equivalent to the standard definition (see example below). However, there may be numerical artifacts in the result that are small in size.
 • Suppose Matrices A and B, respectively, have dimensions (a,b) and (c,d), and define p=a+c-1 and q=b+d-1. Then, the cross-correlation Matrix of dimension (p,q) is given by the formula

${C}_{i,j}={\sum }_{r=1}^{a}\phantom{\rule[-0.0ex]{0.4em}{0.0ex}}{\sum }_{s=1}^{b}{A}_{r,s}\stackrel{&conjugate0;}{{B}_{r-i+c,s-j+d}}$

 where A is assumed to be zero outside of the range of indices (1..a,1..b), and B is assumed to be zero outside of the range of indices (1..c,1..d).
 • FFTs are used when p and q, where p and q are as in the definition above, are powers of 2 and no less than 4 in value. Consequently, if p and q are both a little smaller than a power of 2, the user may want to pad the original Matrices with zeros. See below for an example.
 • If the container=C option is provided, then the results are stored in C. The container must have dimensions (p,q), and datatype complex[8].
 • If A is empty with no elements, a and b are both taken to be zero. Similarly, c and d are both taken to be zero if B has no elements. If both A and B are empty, an error will be thrown because the cross-correlation would have no meaning. However, if, say, B is empty but A is not, then the cross-correlation would be a zero Matrix of size (a-1,b-1).
 • As the underlying implementation of the SignalProcessing package is a module, it is also possible to use the form SignalProcessing:-CrossCorrelation2D to access the command from the package. For more information, see Module Members.

Examples

 > with( SignalProcessing ):

Simple Examples

Example 1

 > A := Vector[row]( [1,2] );
 $\left[\begin{array}{rr}1& 2\end{array}\right]$ (1)
 > B := Matrix( [ [3,-4], [5*I,6] ] );
 $\left[\begin{array}{cc}3& -4\\ 5{}I& 6\end{array}\right]$ (2)
 > C := CrossCorrelation2D( A, B );
 $\left[\begin{array}{ccc}5.999999999999999-1.0877919644084146{}{10}^{-15}{}I& 12.0-5.0{}I& 0.0-9.999999999999998{}I\\ -3.999999999999999-3.6259732146947156{}{10}^{-16}{}I& -4.999999999999999+0.0{}I& 5.999999999999999+0.0{}I\end{array}\right]$ (3)

Example 2

 > A := Matrix( [ [1,2], [3,4] ] );
 $\left[\begin{array}{rr}1& 2\\ 3& 4\end{array}\right]$ (4)
 > B := Matrix( [ [1,2,3], [4,5,6], [7,8,9] ] );
 $\left[\begin{array}{rrr}1& 2& 3\\ 4& 5& 6\\ 7& 8& 9\end{array}\right]$ (5)
 > C1 := CrossCorrelation2D( A, B );
 $\left[\begin{array}{cccc}9.0+0.0{}I& 26.0+0.0{}I& 23.0+0.0{}I& 14.0+0.0{}I\\ 33.0+0.0{}I& 77.0+0.0{}I& 67.0+0.0{}I& 36.0+0.0{}I\\ 21.0+0.0{}I& 47.0+0.0{}I& 37.0+0.0{}I& 18.0+0.0{}I\\ 9.0+0.0{}I& 18.0+0.0{}I& 11.0+0.0{}I& 4.0+0.0{}I\end{array}\right]$ (6)
 > C2 := CrossCorrelation2D( B, A );
 $\left[\begin{array}{cccc}4.0+0.0{}I& 11.0+0.0{}I& 18.0+0.0{}I& 9.0+0.0{}I\\ 18.0+0.0{}I& 37.0+0.0{}I& 47.0+0.0{}I& 21.0+0.0{}I\\ 36.0+0.0{}I& 67.0+0.0{}I& 77.0+0.0{}I& 33.0+0.0{}I\\ 14.0+0.0{}I& 23.0+0.0{}I& 26.0+0.0{}I& 9.0+0.0{}I\end{array}\right]$ (7)

Example 3

 > A := Vector[row]( [I,5,3-I] );
 $\left[\begin{array}{ccc}I& 5& 3-I\end{array}\right]$ (8)
 > B := Vector[row]( [6-0.5*I,7] );
 $\left[\begin{array}{cc}6.-0.5{}I& 7\end{array}\right]$ (9)
 > C := CrossCorrelation2D( A, B );
 $\left[\begin{array}{cccc}0.0+7.0{}I& 34.5+6.0{}I& 51.0-4.5{}I& 18.5-4.5{}I\end{array}\right]$ (10)

Example 4

 > A := Matrix( [ [1,2+7*I], [3-I,4] ] );
 $\left[\begin{array}{cc}1& 2+7{}I\\ 3-I& 4\end{array}\right]$ (11)
 > B := Array( [ [5,6], [7,8+3*I] ] );
 $\left[\begin{array}{cc}5& 6\\ 7& 8+3{}I\end{array}\right]$ (12)
 > C := Matrix( 3, 3, datatype = complex[8] );
 $\left[\begin{array}{ccc}0.0+0.0{}I& 0.0+0.0{}I& 0.0+0.0{}I\\ 0.0+0.0{}I& 0.0+0.0{}I& 0.0+0.0{}I\\ 0.0+0.0{}I& 0.0+0.0{}I& 0.0+0.0{}I\end{array}\right]$ (13)
 > CrossCorrelation2D( A, B, container = C ):
 > 'C' = C;
 ${C}{=}\left(\left[\begin{array}{ccc}8.000000000000004-2.999999999999995{}I& 43.999999999999986+50.0{}I& 14.000000000000004+49.0{}I\\ 27.0-16.999999999999993{}I& 69.99999999999999+23.0{}I& 38.0+34.999999999999986{}I\\ 18.0-5.999999999999995{}I& 38.999999999999986-4.999999999999997{}I& 20.0+7.105427357601002{}{10}^{-15}{}I\end{array}\right]\right)$ (14)

Direct Computation

 • To compute the cross-correlation directly from the definition, first take:
 > X := Matrix( 3, 2, [ [2,-2+I], [-I,2-I], [0,2-2*I] ] );
 $\left[\begin{array}{cc}2& -2+I\\ -I& 2-I\\ 0& 2-2{}I\end{array}\right]$ (15)
 > Y := Matrix( 4, 3, [ [1+I,2+2*I,-1-I], [-2-I,1-I,-2*I], [2*I,1+2*I,1+I], [-1,-I,2-I] ] );
 $\left[\begin{array}{ccc}1+I& 2+2{}I& -1-I\\ -2-I& 1-I& -2{}I\\ 2{}I& 1+2{}I& 1+I\\ -1& -I& 2-I\end{array}\right]$ (16)
 > ( a, b ) := upperbound( X );
 ${a}{,}{b}{≔}{3}{,}{2}$ (17)
 > ( c, d ) := upperbound( Y );
 ${c}{,}{d}{≔}{4}{,}{3}$ (18)
 > p := a + c - 1;
 ${p}{≔}{6}$ (19)
 > q := b + d - 1;
 ${q}{≔}{4}$ (20)
 • We will need to assume X and Y are both 0 when their indices are out of bounds, and take the complex conjugate of the Y values:
 > XX := (i,j) -> if( i < 1 or i > a or j < 1 or j > b, 0, X[i,j] ):
 > YY := (i,j) -> if( i < 1 or i > c or j < 1 or j > d, 0, conjugate( Y[i,j] ) ):
 • The cross-correlation can now be found directly:
 > C1 := Matrix( p, q, proc(i,j) local r, s; add( add( XX(r,s) * YY(r-i+c,s-j+d), s = 1 .. b ), r = 1 .. a ); end proc ):
 • This corresponds to that found with the command:
 > C2 := CrossCorrelation2D( X, Y ):
 > err := max( abs( C1 - C2 ) );
 ${\mathrm{err}}{≔}{3.55271367880050}{}{{10}}^{{-15}}$ (21)

Computation with Fast Fourier Transforms (FFTs)

 • The cross-correlation is computed with FFTs when the dimensions of the final Matrix, namely p and q above, are both powers of 2 and no smaller than 4 in size, and DFTs in the remainder of cases. In the following example, we will compute a cross-correlation both directly (which will use DFTs) and by padding one of the input Matrices (which will use FFTs). The indices 1 and 2 will be used to distinguish the original and padded versions of quantities.
 • First, define the original Matrices:
 > unassign( 'a', 'b', 'c', 'd', 'p', 'q', 'A', 'B', 'C' );
 > a[1] := 60:
 > b[1] := 62:
 > A[1] := LinearAlgebra:-RandomMatrix( a[1], b[1], datatype = complex[8] ):
 > c := 65:
 > d := 60:
 > B := LinearAlgebra:-RandomMatrix( c, d, datatype = complex[8] ):
 • Second, compute the cross-correlation with no padding:
 > C[1] := CrossCorrelation2D( A[1], B ):
 • Now, the dimensions are just a little smaller than powers of 2:
 > p[1] := a[1] + c - 1;
 ${{p}}_{{1}}{≔}{124}$ (22)
 > q[1] := b[1] + d - 1;
 ${{q}}_{{1}}{≔}{121}$ (23)
 • So, we will pad A[1] (padding B would also work) so that the computed cross-correlation Matrix has dimensions that are powers of 2:
 > a[2] := a[1] + 2^(-ilog2(1/max(4,p[1]))) - p[1];
 ${{a}}_{{2}}{≔}{64}$ (24)
 > b[2] := b[1] + 2^(-ilog2(1/max(4,q[1]))) - q[1];
 ${{b}}_{{2}}{≔}{69}$ (25)
 > p[2] := a[2] + c - 1;
 ${{p}}_{{2}}{≔}{128}$ (26)
 > q[2] := b[2] + d - 1;
 ${{q}}_{{2}}{≔}{128}$ (27)
 > A[2] := ArrayTools:-Alias( A[1] ):
 > A[2]( a[1]+1 .. a[2], .. ) := 0:
 > A[2]( .., b[1]+1 .. b[2] ) := 0:
 • Finally, find the cross-correlation using FFTs, and extract the appropriate submatrix:
 > C[2] := CrossCorrelation2D( A[2], B )[ ..p[1], ..q[1] ]:
 • Both results are the same:
 > err := max( abs( C[1] - C[2] ) );
 ${\mathrm{err}}{≔}{8.39987169582269}{}{{10}}^{{-10}}$ (28)

Comparison with 1-D CrossCorrelation

 • The CrossCorrelation2D command can be reduced to the 1-D version with a couple of minor modifications. Take, for example, the following Vectors:
 > A := Vector[row]( [ 5, -2, 3, 4, 1 ] );
 $\left[\begin{array}{rrrrr}5& -2& 3& 4& 1\end{array}\right]$ (29)
 > B := Vector[row]( [ 7, 0, 2 ] );
 $\left[\begin{array}{rrr}7& 0& 2\end{array}\right]$ (30)
 > b := numelems( B );
 ${b}{≔}{3}$ (31)
 • The 2-D cross-correlation returns the following:
 > C__2D := CrossCorrelation2D( A, B );
 $\left[\begin{array}{ccccccc}9.999999999999998-1.6784994417006302{}{10}^{-16}{}I& -4.000000000000003+0.0{}I& 41.0+0.0{}I& -5.999999999999999+0.0{}I& 22.999999999999993+0.0{}I& 27.999999999999996+0.0{}I& 7.000000000000003+0.0{}I\end{array}\right]$ (32)
 • The 1-D cross-correlation, however, gives:
 > CrossCorrelation( A, B );
 $\left[\begin{array}{ccccccc}41.0& -4.0& 10.0& 0.0& 0.0& 0.0& 0.0\end{array}\right]$ (33)
 • If we switch the arguments for CrossCorrelation, and choose the appropriate value of lowerlag, the two cross-correlations agree:
 > C__1D := CrossCorrelation( B, A, 1-b );
 $\left[\begin{array}{ccccccc}10.0& -4.0& 41.0& -6.0& 23.0& 28.0& 7.0\end{array}\right]$ (34)

Application to Fingerprint Matching

 • A classic application of two-dimensional cross-correlation is the lining up of two or more images, for example, images of fingerprint fragments. Consider the following image:
 > with( ImageTools ):
 > imgfile := cat( kernelopts( mapledir ), kernelopts( dirsep ), data, kernelopts( dirsep ), "images", kernelopts( dirsep ), "fingerprint.jpg" ):
 > img1 := Matrix( Read( imgfile ) ):
 > Preview( img1 );
 • First, determine the dimensions:
 > ( a, b ) := upperbound( img1 );
 ${a}{,}{b}{≔}{240}{,}{256}$ (35)
 • We will also use the mean later:
 > mu := add( img1 ) / numelems( img1 ):
 • Second, let's choose some bounds for a fragment of the fingerprint, extract the subimage (fingerprint fragment), and then add some noise (to make our task a little more realistic):
 > cL := 85:
 > cU := 205:
 > dL := 120:
 > dU := 180:
 > img2 := img1[ cL..cU, dL..dU ] + 0.5 * mu * ArrayTools:-RandomArray( cU-cL+1, dU-dL+1, distribution = normal ):
 > Preview( img2 );
 • Suppose now that we do not know whether the fragment is part of the larger image, and wish to determine if the fragment is a "match". To do this, compute the cross-correlation of the tempered data (formed by subtracting the mean from all elements):
 > mu := add( img1 ) / numelems( img1 ):
 > C := Matrix( Re( CrossCorrelation2D( img1 -~ mu, img2 -~ mu ) ), datatype = sfloat ):
 • The maximum correlation occurs here, which, as expected, coincides with the bottom-right corner of the fragment:
 > ( alpha, beta ) := max[index]( C );
 ${\mathrm{\alpha }}{,}{\mathrm{\beta }}{≔}{205}{,}{180}$ (36)
 • Visually:
 > matplot := plots:-matrixplot( C ):
 > plots:-display( matplot, title = "Cross-Correlation - Full" );
 > r := 25:
 > box := [ alpha-r .. alpha+r, beta-r .. beta+r, min(C) .. max(C) ]:
 > plots:-display( matplot, view = box, title = "Cross-Correlation - Zoomed" );
 >

Compatibility

 • The SignalProcessing[CrossCorrelation2D] command was introduced in Maple 2020.
 • For more information on Maple 2020 changes, see Updates in Maple 2020.