John Hunter writes:
"Todd" == Todd Miller <jmiller@stsci.edu> writes:
>> I have a series of luminance images that I want to do some >> correlation analyses on. Each image is an MxN frame of a >> movie. I want to compute the correlation between a given pixel >> i,j, and every other pixel in the image over each frame. That >> is, if I assume xij is a numFrames length time series, and xkl >> is another numFrames length time series (the pixel intensities >> at two points in the image over time), I want to compute the
>> corrcoeff(xij, xkl) for every kl with ij fixed.
>> I know I could do this by looping over the pixels of the image, >> but I'm hoping for something a bit faster.
Todd> For numarray try numarray.convolve.correlate2d and set Todd> fft=1.
Something that Todd and I have talked about is making general functions (like fft, convolve, correlate) broadcastable to other dimensions. Not much has been done to do that but I think a general mechanism could be developed to do just that. In that way one could do a 1-d correlation over a 2 or 3-d array without looping over the extra dimensions (all that would be done implicitly in C) with the option of specifying which dimension the function should be applied to while looping over all the other dimensions. In the meantime, it would seemt that one possible way of dealing with this is to simply recast your array as a 1-d array to use the 1-d correlation. This means some memory copying and explicit padding (how do you want the correlation to handle points beyond the ends?). Supposing you want it to wraparound (in effect a circular correlation) you could do something like this (untested and probably has some mistakes :-) if imagecube is a T x M x N array where there are T time samples. data = concatenate((imagecube, imagecube, imagecube)) # pad the time series on both ends with identical copies reference = imagecube[:, i, j] flatdata = ravel(transpose(data)) flatresult = correlate(flatdata, reference) result = transpose(reshape(flatresult, (N,M,T)))[T:2*T,:,:] This is admittedly a bit clumsy, but should be much, much faster than iterating over all pixels. It would be much nicer just to have result = correlate(imagecube, imagecube[:,i,j], axis=0) which the broadcasting facility would permit. Perry