strange behavior convolving via fft
at least I think this is strange behavior. When convolving an image with a large kernel, its know that its faster to perform the operation as multiplication in the frequency domain. The below code example shows that the results of my 2d filtering are shifted from the expected value a distance 1/2 the width of the filter in both the x and y directions. Can anyone explain why this occurs? I have been able to find the answer in any of my image processing books. The code sample below is an artificial image of size (100, 100) full of zeros, the center of the image is populated by a (10, 10) square of 1's. The filter kernel is also a (10,10) square of 1's. The expected result of the convolution would therefore be a peak at location (50,50) in the image. Instead, I get (54, 54). The same shifting occurs regardless of the image and filter (assuming the filter is symetric, so flipping isnt necessary). I came across this behavior when filtering actual images, so this is not a byproduct of this example. The same effect also occurs using the full FFT as opposed to RFFT. I have links to the images produced by this process below. Thanks for any insight anyone can give! Chris In [12]: a = np.zeros((100,100)) In [13]: a[45:55,45:55] = 1 In [15]: k = np.ones((10,10)) In [16]: afft = np.fft.rfft2(a, s=(256,256)) In [19]: kfft = np.fft.rfft2(k, s=(256,256)) In [21]: result = np.fft.irfft2(afft*kfft).real[0:100,0:100] In [23]: result.argmax() Out[23]: 5454 www.therealstevencolbert.com/dump/a.png www.therealstevencolbert.com/dump/afft.png www.therealstevencolbert.com/dump/kfft.png www.therealstevencolbert.com/dump/result.png
On Mon, May 11, 2009 at 9:40 AM, Chris Colbert <sccolbert@gmail.com> wrote:
at least I think this is strange behavior.
When convolving an image with a large kernel, its know that its faster to perform the operation as multiplication in the frequency domain. The below code example shows that the results of my 2d filtering are shifted from the expected value a distance 1/2 the width of the filter in both the x and y directions. Can anyone explain why this occurs? I have been able to find the answer in any of my image processing books.
The code sample below is an artificial image of size (100, 100) full of zeros, the center of the image is populated by a (10, 10) square of 1's. The filter kernel is also a (10,10) square of 1's. The expected result of the convolution would therefore be a peak at location (50,50) in the image. Instead, I get (54, 54). The same shifting occurs regardless of the image and filter (assuming the filter is symetric, so flipping isnt necessary).
Your kernel is offset and the result is expected. The kernel needs to be centered on the origin, aliasing will then put parts of it in all four corners of the array *before* you transform it. If you want to keep it simple you can phase shift the transform instead. Chuck
Ok, that makes sense. Thanks Chuck. On Mon, May 11, 2009 at 2:41 PM, Charles R Harris <charlesr.harris@gmail.com
wrote:
On Mon, May 11, 2009 at 9:40 AM, Chris Colbert <sccolbert@gmail.com>wrote:
at least I think this is strange behavior.
When convolving an image with a large kernel, its know that its faster to perform the operation as multiplication in the frequency domain. The below code example shows that the results of my 2d filtering are shifted from the expected value a distance 1/2 the width of the filter in both the x and y directions. Can anyone explain why this occurs? I have been able to find the answer in any of my image processing books.
The code sample below is an artificial image of size (100, 100) full of zeros, the center of the image is populated by a (10, 10) square of 1's. The filter kernel is also a (10,10) square of 1's. The expected result of the convolution would therefore be a peak at location (50,50) in the image. Instead, I get (54, 54). The same shifting occurs regardless of the image and filter (assuming the filter is symetric, so flipping isnt necessary).
Your kernel is offset and the result is expected. The kernel needs to be centered on the origin, aliasing will then put parts of it in all four corners of the array *before* you transform it. If you want to keep it simple you can phase shift the transform instead.
Chuck
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Hi Chris 2009/5/11 Chris Colbert <sccolbert@gmail.com>:
When convolving an image with a large kernel, its know that its faster to perform the operation as multiplication in the frequency domain. The below code example shows that the results of my 2d filtering are shifted from the expected value a distance 1/2 the width of the filter in both the x and y directions. Can anyone explain why this occurs? I have been able to find the answer in any of my image processing books.
Just as a reminder, when doing this kind of filtering always pad correctly. Scipy does this in scipy.signal.fftconvolve I've also got some filtering implemented in http://mentat.za.net/cgi-bin/hgwebdir.cgi/filter/file/e97c0a6dd0ea/lpi_filte... Regards Stéfan
Stefan, Did I pad my example incorrectly? Both images were upped to the larger nearest power of 2 (256)... Does the scipy implementation do this differently? I thought that since FFTW support has been dropped, that scipy and numpy use the same routines... Thanks! Chris 2009/5/11 Stéfan van der Walt <stefan@sun.ac.za>
Hi Chris
When convolving an image with a large kernel, its know that its faster to perform the operation as multiplication in the frequency domain. The below code example shows that the results of my 2d filtering are shifted from
expected value a distance 1/2 the width of the filter in both the x and y directions. Can anyone explain why this occurs? I have been able to find
2009/5/11 Chris Colbert <sccolbert@gmail.com>: the the
answer in any of my image processing books.
Just as a reminder, when doing this kind of filtering always pad correctly. Scipy does this in
scipy.signal.fftconvolve
I've also got some filtering implemented in
http://mentat.za.net/cgi-bin/hgwebdir.cgi/filter/file/e97c0a6dd0ea/lpi_filte...
Regards Stéfan _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Hi Chris, If you have MxN and PxQ signals, you must pad them to shape M+P-1 x N+Q-1, in order to prevent circular convolution (i.e. values on the one end sliding back in at the other). Regards Stéfan 2009/5/11 Chris Colbert <sccolbert@gmail.com>:
Stefan,
Did I pad my example incorrectly? Both images were upped to the larger nearest power of 2 (256)...
Does the scipy implementation do this differently? I thought that since FFTW support has been dropped, that scipy and numpy use the same routines...
Thanks!
Chris
2009/5/11 Chris Colbert <sccolbert@gmail.com>:
Does the scipy implementation do this differently? I thought that since FFTW support has been dropped, that scipy and numpy use the same routines...
Just to be clear, I was referring to scipy.signal.fftconvolve, not scipy's FFT (which is the same as NumPy's). Regards Stéfan
Thanks Stefan. 2009/5/11 Stéfan van der Walt <stefan@sun.ac.za>
2009/5/11 Chris Colbert <sccolbert@gmail.com>:
Does the scipy implementation do this differently? I thought that since FFTW support has been dropped, that scipy and numpy use the same routines...
Just to be clear, I was referring to scipy.signal.fftconvolve, not scipy's FFT (which is the same as NumPy's).
Regards Stéfan _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
participants (3)
-
Charles R Harris -
Chris Colbert -
Stéfan van der Walt