I am looking for efficient ways to code neighbourhood functions. For example a neighbourhod add for an element in an array will simply be the sum of the neighbours: 1 0 2 3 x 3 , then x becomes 7 (first order neighbour), 11 (2nd order) etc. 1 1 0 I would be interested in efficient ways of doing this for a whole array, something like a_nsum = neighbour_sum(a, order=1), where each element in a_nsum is the sum of the corresponding element in a. There must be some work done on neighbourhood functions for arrays, so I would be grateful for some pointers. Thanks, Robert Denham ************************************************************************ The information in this e-mail together with any attachments is intended only for the person or entity to which it is addressed and may contain confidential and/or privileged material. Any form of review, disclosure, modification, distribution and/or publication of this e-mail message is prohibited. If you have received this message in error, you are asked to inform the sender as quickly as possible and delete this message and any copies of this message from your computer and/or your computer system network. ************************************************************************
Hi Robert, Isn't this just a convolution? For your first order case you are convolving the array with a kernel that looks like: 0 1 0 1 0 1 0 1 0 and for second order: 1 1 1 1 0 1 1 1 1 If this is what you are looking for, there is a convolve function built in to Numeric but it is for rank-1 arrays only. You can use 2D FFTs for a 2D case -- although the efficiency won't be the greatest unless you can use the FFT of the kernel array and/or the data array over and over again (since in general a convolution by FFTs takes 3 FFTs -- 1 for the data, 1 for the kernel, and 1 inverse one after you have multiplied the first two together). It would work well for higher order cases...(but beware of the wrapping that goes on accross the boundaries!) For lower order cases or when you can't re-use the FFTs, you'll probably want a brute force technique -- which I'll leave for someone else... Scott Robert.Denham@dnr.qld.gov.au wrote:
-- Scott M. Ransom Address: Harvard-Smithsonian CfA Phone: (617) 495-4142 60 Garden St. MS 10 email: ransom@cfa.harvard.edu Cambridge, MA 02138 GPG Fingerprint: 06A9 9553 78BE 16DB 407B FFCA 9BFA B6FF FFD3 2989
As Scott Ransom already mentioned, these "neighborhood functions" are usually referred to as convolutions and are widely used in signal/image processing. For large convolution kernels, the most efficient implementation is to use Fourier methods, since a convolution in the spatial domain maps to a multiplication in the Fourier domain. However for small kernels this is inefficient because the time spent doing the forward and inverse FFT's dwarfs the time that it would take to just do the convolution. There is a convolve function built into Numeric but it only is implemented for 1-d arrays. It would be nice if this were generalized... when somebody gets the time. In the meanwhile - here are two more comments which may help. If your kernel is separable (i.e. a rank-one matrix, or equivalently, the outer product of a column and a row vector) then the 2-d convolution is equivalent to doing 2 1-d convolutions in sequence. For your "first order neighborhood function" the kernel is 0 1 0 1 0 1 0 1 0 which is not separable. But for the "second order" case, the kernel is 1 1 1 1 0 1 1 1 1 which is *almost* separable if you made that middle 0 into a 1. But if you were to convolve with the kernel 1 1 1 1 1 1 1 1 1 then subtract off the value of the original array, then you'd have what you were looking for. And convolving with this kernel is essentially the same as convolving with the 1-d kernel [1 1 1], then transposing, then convolving with [1 1 1] again, and transposing a second time. This scales up to larger separable kernels. I'm not sure how efficient this winds up being - transposing large arrays can be slow, if the arrays are too large to sit in physical memory - the access pattern of a transpose is hell on virtual memory. Finally, I would suggest that you look at the file Demo/life.py in the Numeric Python distribution - in particular the functions shift_up, shift_down, shift_left, and shift_right. Using these you could write: def first_order(arr): return shift_left(arr) + shift_right(arr) + shift_up(arr) + shift_down(arr) which is nice and terse and understandable, but unfortunately not very memory-efficient. If speed and memory usage are both critical, your best bet is probably to write the functions in C and use SWIG to create Python bindings. Charles G. Waldman Object Environments cgw@objenv.com
Hi Robert, Isn't this just a convolution? For your first order case you are convolving the array with a kernel that looks like: 0 1 0 1 0 1 0 1 0 and for second order: 1 1 1 1 0 1 1 1 1 If this is what you are looking for, there is a convolve function built in to Numeric but it is for rank-1 arrays only. You can use 2D FFTs for a 2D case -- although the efficiency won't be the greatest unless you can use the FFT of the kernel array and/or the data array over and over again (since in general a convolution by FFTs takes 3 FFTs -- 1 for the data, 1 for the kernel, and 1 inverse one after you have multiplied the first two together). It would work well for higher order cases...(but beware of the wrapping that goes on accross the boundaries!) For lower order cases or when you can't re-use the FFTs, you'll probably want a brute force technique -- which I'll leave for someone else... Scott Robert.Denham@dnr.qld.gov.au wrote:
-- Scott M. Ransom Address: Harvard-Smithsonian CfA Phone: (617) 495-4142 60 Garden St. MS 10 email: ransom@cfa.harvard.edu Cambridge, MA 02138 GPG Fingerprint: 06A9 9553 78BE 16DB 407B FFCA 9BFA B6FF FFD3 2989
As Scott Ransom already mentioned, these "neighborhood functions" are usually referred to as convolutions and are widely used in signal/image processing. For large convolution kernels, the most efficient implementation is to use Fourier methods, since a convolution in the spatial domain maps to a multiplication in the Fourier domain. However for small kernels this is inefficient because the time spent doing the forward and inverse FFT's dwarfs the time that it would take to just do the convolution. There is a convolve function built into Numeric but it only is implemented for 1-d arrays. It would be nice if this were generalized... when somebody gets the time. In the meanwhile - here are two more comments which may help. If your kernel is separable (i.e. a rank-one matrix, or equivalently, the outer product of a column and a row vector) then the 2-d convolution is equivalent to doing 2 1-d convolutions in sequence. For your "first order neighborhood function" the kernel is 0 1 0 1 0 1 0 1 0 which is not separable. But for the "second order" case, the kernel is 1 1 1 1 0 1 1 1 1 which is *almost* separable if you made that middle 0 into a 1. But if you were to convolve with the kernel 1 1 1 1 1 1 1 1 1 then subtract off the value of the original array, then you'd have what you were looking for. And convolving with this kernel is essentially the same as convolving with the 1-d kernel [1 1 1], then transposing, then convolving with [1 1 1] again, and transposing a second time. This scales up to larger separable kernels. I'm not sure how efficient this winds up being - transposing large arrays can be slow, if the arrays are too large to sit in physical memory - the access pattern of a transpose is hell on virtual memory. Finally, I would suggest that you look at the file Demo/life.py in the Numeric Python distribution - in particular the functions shift_up, shift_down, shift_left, and shift_right. Using these you could write: def first_order(arr): return shift_left(arr) + shift_right(arr) + shift_up(arr) + shift_down(arr) which is nice and terse and understandable, but unfortunately not very memory-efficient. If speed and memory usage are both critical, your best bet is probably to write the functions in C and use SWIG to create Python bindings. Charles G. Waldman Object Environments cgw@objenv.com
participants (3)
-
Charles G Waldman -
Robert.Denham＠dnr.qld.gov.au -
Scott Ransom