[SciPy-User] moving window (2D) correlation coefficient

[Vincent Schut]

On 02/10/2010 05:42 PM, Zachary Pincus wrote:
> I bet that you could construct an array with shape (x,y,5,5), where
> array[i,j,:,:] would give the 5x5 window around (i,j), as some sort of
> mind-bending view on an array of shape (x,y), using a positive offset
> and some dimensions having negative strides.

I tried that, but then less clever (just reshaping and transposing); it 
made a copy, and I ended up with an array that took 5*5 times as much 
space as my initial array... And it took a lot of time.
But I don't see how this reshape would gain me much speed?

  Then you could compute
> the correlation coefficient between the two arrays directly. Maybe?

Ah, like that. I thought numpy.corrcoef(x,y) only worked on flat x and 
flat y. There is no axis keyword... Is there another function that would 
work on nd arrays?

>
> Probably cython would be more maintainable...

Well, speed is more important than readability this time :-) It's 
terabytes of data I'll need to push through this function...
What I understood till now (just starting to look at cython) is that in 
this kind of thing, the repeated array sub-indexing (and its bounds 
checking etc.) is the main bottleneck. I presume the corrcoef function 
is pretty much optimized? So I thought if I could do the loop in cython 
and implement some clever code for the window indexing (only replace the 
new elements), I'd relatively simple get the major speed gain. Does that 
sound correct?
Only still have to find out now how to call the corrcoef function from 
cython...

Thanks for your thoughts on this,
Vincent.
>
> Zach
>
>
> On Feb 10, 2010, at 10:45 AM, Vincent Schut wrote:
>
>> Hi,
>>
>> I need to calculate the correlation coefficient of a (simultaneously)
>> moving window over 2 images, such that the resulting pixel x,y (center
>> of the window) is the corrcoef((a 5x5 window around x,y for image
>> A), (a
>> 5x5 window around x,y for image B)).
>> Currently, I just loop over y, x, and then call corrcoef for each x,y
>> window. Would there be a better way, other than converting the loop to
>> cython?
>>
>>
>> For clarity (or not), the relevant part from my code:
>>
>>
>> for y in range(d500shape[2]):
>>      for x in range(d500shape[3]):
>>          if valid500[d,y,x]:
>>              window = spectral500Bordered[d,:,y:y+5, x:x
>> +5].reshape(7, -1)
>>              for b in range(5):
>>                  nonzeroMask = (window[0]>  0)
>>                  b01corr[0,b,y,x] =
>> numpy.corrcoef(window[0][nonzeroMask], window[b+2][nonzeroMask])[0,1]
>>                  b01corr[1,b,y,x] =
>> numpy.corrcoef(window[1][nonzeroMask], window[b+2][nonzeroMask])[0,1]
>>
>>
>> forget the 'if valid500' and 'nonzeroMask', those are to prevent
>> calculating pixels that don't need to be calculated, and to feed only
>> valid pixels from the window into corrcoef
>> spectral500Bordered is essentially a [d dates, 7 images, ysize, xsize]
>> array. I work per date (d), then calculate the corrcoef for images[0]
>> versus images[2:], and for images[1] versus images[2:]
>>
>> Thanks,
>> Vincent.
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user