[SciPy-User] moving window (2D) correlation coefficient
josef.pktd@gmai...
josef.pktd@gmai...
Wed Feb 10 13:48:23 CST 2010
On Wed, Feb 10, 2010 at 2:36 PM, Vincent Schut <schut@sarvision.nl> wrote:
> On 02/10/2010 05:53 PM, josef.pktd@gmail.com wrote:
>> On Wed, Feb 10, 2010 at 11:42 AM, Zachary Pincus
>> <zachary.pincus@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@scipy.org
>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User@scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
More information about the SciPy-User
mailing list