-2
$\begingroup$

I've defined a correlation function which does the following:

  • Take two data sets $d_1$ and $d_2$ which are the same size and shape
  • Pad $d_2$ with zeroes (to avoid non-positive indices in the next step)
  • Crop $d_2$ to a region the size and shape of the original image but displaced from the centre by $r_n$ rows and $c_n$ columns
  • Multiply $d_1\,\!^\ast$ by the displaced $d_2$
  • Normalise the result by overlapping area (to avoid the possibility that a large displacement will result in a small correlation simply because there are not enough pixels rather than because the correlation is actually low)

Here's the code ($rows$ is the number of rows in each data set, and $cols$ is the number of columns in each data set):

corr[d1_,d2_,rn_,cn_]:=Total[
                        Total[
                         Conjugate[d1]ArrayPad[
                          d2,{{Length[d2]},{Length[Transpose[d2]]}}
                         ][[Length[d1]+1+rn;;2Length[d1]+rn,
                            Length[Transpose[d1]]+1+cn;;2Length[Transpose[d1]]+cn]]
                        ]
                       ]/((rows-rn)(cols-cn));

(Before anyone says anything, yes, I do know Mathematica already has a built-in correlation function. I made my own correlation function because I didn't want to spend a long time checking that Mathematica's function did exactly what I wanted and possibly troubleshooting it. I create a lot of zeroes, sure, but in the end I'm just multiplying $rows\times cols$ pairs of numbers and adding the results together.)

Given two data sets, I then do this for different values of $r_n$ and $c_n$ and create a correlation map:

corrdata=ParallelTable[corr[data1,data2,rn,cn],{rn,-100,100},{cn,-150,150}]

My problem is that I have big data sets and this takes a long time. I want to run this correlation function on 10,000 data sets. If running it on data sets with 80,000 pixels and making correlation maps with 60,501 pixels takes 15 minutes on my computer, running it 10,000 times is going to take 104 days.

Is there a way to make this (a lot) quicker?

Thanks in advance.

$\endgroup$
  • $\begingroup$ For a start, Length[Transpose[d2] should be replaced by Dimensions[d1][[2]]2. Or would you turn an elefant onto its back just for counting its legs? $\endgroup$ – Henrik Schumacher May 8 '18 at 11:54
  • $\begingroup$ Don't use ArrayPad; instead extract the correct rows and colums of d1 and d2 with Part. Having extracted them into, say a1 and a2, use Tr[Conjugate[a1].a2] or Flatten[Conjugate[a1]].Flatten[a2] instead multiplying and Total. $\endgroup$ – Henrik Schumacher May 8 '18 at 11:58
  • 3
    $\begingroup$ Why not use the built in function ListCorrelate? Since you are worried about efficiency, there's a good chance it will be faster. $\endgroup$ – bill s May 8 '18 at 12:04
  • $\begingroup$ Edit: Should be Tr[ConjugateTranspose[a1].a2], but Flattening is problably faster. $\endgroup$ – Henrik Schumacher May 8 '18 at 12:07
  • 3
    $\begingroup$ I think just ListCorrelate[d1, d2, {-1, 1}, 0] would be orders of magnitude faster, although you will have to also determine the normalization, which can be determined by replacing d1 and d2 with matrices of 1s. $\endgroup$ – Carl Woll May 8 '18 at 15:12
4
$\begingroup$

I don't know why you insist on not using ListCorrelate. Here is a version of your code that uses ListCorrelate:

listC[d1_,d2_]:=Module[{r, c},
    {r,c}=Dimensions[d1];
    Divide[
        ListCorrelate[Conjugate[d1],d2,{-1,1},0],
        N@Outer[Times,Reverse@Range[2r-1],Reverse@Range[2c-1]]
    ]
]

(by the way, I think your computation of the denominator is flawed, but I use the same denominator so that the outputs can be compared). And, here is a comparison with your code in the OP:

d1 = RandomComplex[1+I, {20,30}];
d2 = RandomComplex[1+I, {20,30}];

r1 = listC[d1, d2]; //AbsoluteTiming
r2 = Table[corr[d1,d2,rn,cn], {rn,-19,19}, {cn,-29,29}]; //AbsoluteTiming

Block[{Internal`$EqualTolerance=4}, r1 == r2]

{0.000654, Null}

{1.13723, Null}

True

(I needed to lower the tolerance due to numeric fuzz). As you can see, ListCorrelate is already over 1000 times faster. Here is an example using larger matrices:

d1 = RandomComplex[1+I, {1000,1000}];
d2 = RandomComplex[1+I, {1000,1000}];

listC[d1, d2]; //AbsoluteTiming

{0.704561, Null}

$\endgroup$
  • $\begingroup$ I don't know why people insist on using Mathematica's built-in functions. I clearly stated my reasons not to want to use them right from the beginning, yet I keep receiving suggestions to use them. I appreciate the help, but not the insistence that I do things I've already said I don't want to do. $\endgroup$ – Rain May 8 '18 at 23:09
  • $\begingroup$ My denominator is indeed wrong. I should have used Abs[rn] and Abs[cn] rather than rn and cn. Thanks for pointing it out. $\endgroup$ – Rain May 8 '18 at 23:09
  • $\begingroup$ Thanks for the answer. I'll try ListCorrelate after my current code finishes running. $\endgroup$ – Rain May 8 '18 at 23:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.