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.comwrote:
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
Numpydiscussion mailing list Numpydiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
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/cgibin/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
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/cgibin/hgwebdir.cgi/filter/file/e97c0a6dd0ea/lpi_filte...
Regards Stéfan _______________________________________________ Numpydiscussion mailing list Numpydiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
Hi Chris,
If you have MxN and PxQ signals, you must pad them to shape M+P1 x N+Q1, 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 _______________________________________________ Numpydiscussion mailing list Numpydiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
participants (3)

Charles R Harris

Chris Colbert

Stéfan van der Walt