Helder writes:
>
> On Wednesday, November 6, 2013 4:13:49 PM UTC+1, Mats Löfdahl wrote:
>> I found something surprising (to me) with the convolve() IDL function. There is something strange about how it does its Fourier wrap-around of an image from one side of the array to the other.
>>
>>
>>
>> Here is a simple example. First define a simple image where half is unity and half is zero:
>>
>>
>>
>> sz = 10
>>
>> im = [replicate(1., sz/2), replicate(0., sz/2)] # replicate(1., sz)
>>
>> print, 'Original:'
>>
>> print,im[*,sz/2], format = '(f5.2)'
>>
>>
>>
>> This gives the output:
>>
>>
>>
>> Original:
>>
>> 1.00
>>
>> 1.00
>>
>> 1.00
>>
>> 1.00
>>
>> 1.00
>>
>> 0.00
>>
>> 0.00
>>
>> 0.00
>>
>> 0.00
>>
>> 0.00
>>
>>
>>
>>
>>
>> Then define a point spread function and do the convolution:
>>
>>
>>
>> psf1 = [[1., 1., 1.], [1., 5., 1.], [1., 1., 1.]]
>>
>> psf1 = psf1/total(psf1)
>>
>> imc1 = convolve(im, psf1)
>>
>> print, 'With convolve:'
>>
>> print,imc1[*,sz/2], format = '(f5.2)'
>>
>>
>>
>> The output I get is:
>>
>>
>>
>> With convolve:
>>
>> 0.77
>>
>> 1.00
>>
>> 1.00
>>
>> 1.00
>>
>> 0.77
>>
>> 0.23
>>
>> -0.00
>>
>> 0.00
>>
>> -0.00
>>
>> -0.00
>>
>>
>>
>> See how the wrap-around reduced the 1.00 in the first pixel to 0.75 but the last pixel does not get the corresponding increase?
>>
>>
>>
>> Whereas if I do the equivalent operation explicitly with FFT, I do get the expected 0.23 in the last pixel:
>>
>>
>>
>> psf2 = fltarr(sz, sz)
>>
>> psf2[sz/2-1:sz/2+1, sz/2-1:sz/2+1] = psf1*sz*sz
>>
>> psf2 = shift(psf2, sz/2, sz/2)
>>
>> imc2 = float(fft(fft(im)*fft(psf2), /inv))
>>
>> print, 'With fft:'
>>
>> print,imc2[*,sz/2], format = '(f5.2)'
>>
>>
>>
>> With fft:
>>
>> 0.77
>>
>> 1.00
>>
>> 1.00
>>
>> 1.00
>>
>> 0.77
>>
>> 0.23
>>
>> -0.00
>>
>> -0.00
>>
>> -0.00
>>
>> 0.23
>>
>>
>>
>> I've looked at the code in http://www.astro.washington.edu/docs/idl/cgi-bin/getpro/libr ary21.html?CONVOLVE and as far as I can see (due to various options the code is not entirely straight forward to read), the fft convolution has no reason to do give any different result from what I do explicitly with fft.
>>
>>
>>
>> Does anybody know what is going on?
>
> Dunno,
> but if you try this, you get what you expected:
>
> imc2 = CONVOL_FFT(im, psf1, /NO_PADDING)
> print,imc2[*,sz/2], format = '(f5.2)'
> 0.77
> 1.00
> 1.00
> 1.00
> 0.77
> 0.23
> -0.00
> -0.00
> -0.00
> 0.23
>
> but of course you need at least IDL 8.1.
> With padding you get the "wrong" result.
You get the result you expect if you use the IDL routine CONVOL (instead
of the CONVOLVE you are using, and set the EDGE_WRAP keyword:
psf1 = [[1., 1., 1.], [1., 5., 1.], [1., 1., 1.]]
psf1 = psf1/total(psf1)
imc1 = convol(im, psf1, /edge_wrap)
print, 'With convol:'
print,imc1[*,sz/2], format = '(f5.2)'
END
IDL> .go
% Compiled module: $MAIN$.
With convol:
0.77
1.00
1.00
1.00
0.77
0.23
0.00
0.00
0.00
0.23
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
Sepore ma de ni thue. ("Perhaps thou speakest truth.")
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
