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

[Vincent Schut]

On 02/10/2010 05:53 PM, josef.pktd at gmail.com wrote:
> On Wed, Feb 10, 2010 at 11:42 AM, Zachary Pincus
> <zachary.pincus at yale.edu>  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. Then you could compute
>> the correlation coefficient between the two arrays directly. Maybe?
>>
>> Probably cython would be more maintainable...
>>
>> 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:]
>
> I wrote a moving correlation for time series last november (scipy user
> and preceding discussion on numpy mailing list)
> I don't work with pictures, so I don't know if this can be extended to
> your case. Since signal.correlate or convolve work in all directions
> it might be possible
>
>     def yxcov(self):
>         xys = signal.correlate(self.x*self.y, self.kern,
> mode='same')[self.sslice]
>         return xys/self.n - self.ymean*self.xmean
>
> Josef

I saw that when searching on this topic, but didn't think it would work 
for me as I supposed it was purely 1-dimensional, and I thought that in 
your implementation, though the window moves, the kernel is the same all 
the time? I'm no signal processing pro (alas) so please correct me if 
I'm wrong. I'll try to find the discussion you mentioned tomorrow. Damn 
time zones... ;-)
>
>>>
>>> Thanks,
>>> Vincent.
>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>