[Numpy-discussion] image correlations
perry at stsci.edu
Tue May 25 11:52:01 CDT 2004
John Hunter writes:
> >>>>> "Todd" == Todd Miller <jmiller at 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.
More information about the Numpy-discussion