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

[Vincent Schut]

On 02/10/2010 10:24 PM, David Baddeley wrote:
> I think Josef's nailed it - if you can live with subtracting a moving mean from each pixel [calculated using a correlation/convolution with (1.0/(windowsize*windowsize))*ones((windowsize, windowsize))] rather than taking the mean across the actual window used for the correlation, you'll be home free. Not only do you save yourself the loops, but you're also going to be doing significantly less multiplications than if you calculated the correlation separately for each window.
>
> I suspect that this approach would even be considerably faster than coding your loops and the windowed correlation up in cython.
>
> cheers,
> David
>

Guys, thanks a lot! It works well, and executing times have decreased 
with approx a factor 6. Pity I still don't get to learn cython, but I 
gained some insight. Be asured I'd buy you a beer if we lived a bit 
closer :-)

thanks again,
Vincent.

>
> ----- Original Message ----
> From: "josef.pktd at gmail.com"<josef.pktd at gmail.com>
> To: SciPy Users List<scipy-user at scipy.org>
> Sent: Thu, 11 February, 2010 8:48:23 AM
> Subject: Re: [SciPy-User] moving window (2D) correlation coefficient
>
> On Wed, Feb 10, 2010 at 2:36 PM, Vincent Schut<schut at sarvision.nl>  wrote:
>> 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... ;-)
>
> correlation is just  sum(x*y)/sqrt(sum(x*x))/sqrt(sum(y*y))   but
> subtracting moving mean which makes it slightly tricky.
>
> the moving sum can be done with any nd convolve or correlate (ndimage
> or signal) which goes in all directions
>
> the window is just np.ones((windowsize,windowsize)) for symmetric
> picture, or if nan handling is not necessary than the averaging window
> can be used directly.
>
> I'm pretty sure now it works, with 5 or so intermediate arrays in the
> same size as the original array
>
> Josef
>
>
>
>>>
>>>>>
>>>>> 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
>>>>
>>
>>
>> _______________________________________________
>> 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