Galsim Convolution vs scipy fftconvolve

[jecampagne]

Hi,
I was questioning the Galsim convolution algo and I have made a simple exo using version 2.4.8 (the Galsim version install was done through JAX-Galsim as the direct install of GalSim in Colab fails, see my previous issue)

gal_flux = 1.e5    # total counts on the image
gal_sigma = 2.     # arcsec
psf_sigma = 1.     # arcsec
pixel_scale = 0.2  # arcsec / pixel

nx,ny=128,128

gal = _galsim.Gaussian(flux=gal_flux, sigma=gal_sigma)
psf = _galsim.Gaussian(flux=1., sigma=psf_sigma)
final = _galsim.Convolve([gal, psf])

method =  'auto' # 'real_space' 'auto' 'fft' (the result does not change wrt to the method s I thing FFT is used at the end)
im_final = final.drawImage(nx=nx,ny=ny,scale=pixel_scale, method=method)

from scipy.signal import fftconvolve as sc_fftconvolve
sci_conv2d = sc_fftconvolve(im_gal.array,im_psf.array, mode="same")
#plot

Does one fully understand the difference between the convolution done by Galsim and FFTconvolve ?

Thanks
JE

PS: the effect of the Pixel convolution is 2 order of magnitudes smaller than the difference between Galsim and scipy fftconvolve.

[rmandelb]

To me the difference image looks like a difference in centroiding conventions of some kind. Rather than a purely visual comparison, have you run galsim.hsm.FindAdaptiveMom to quantitatively compare the moments (0th, 1st, and 2nd) in the two cases? That would help confirm whether it's purely a centroiding difference or something else.

As you noted, the other difference could be pixel response. I don't see in that code snippet where you've drawn im_gal and im_psf, but if each includes the pixel response, then pixel response is included in the GalSim final object once, and in the scipy convolution twice. But that difference would be easy to eliminate by drawing one of the objects without including the pixel response.

[jmeyers314]

Seems like some kind of indexing or centroiding difference-in-assumption from scipy.
If you ask what the peak pixel is in each of the inputs via something likenp.unravel_index(np.argmax(array), array.shape), you'll find the peak is (63, 63) for the galaxy image, psf image, and galsim convolution image (as makes sense for these circular monotonically declining profiles). But the scipy result has its peak at (64, 64). I tried rolling the scipy image around by 1 pixel in various directions, but this didn't seem to actually help, just moved the dipole around in different directions or made it stronger so must be a smaller-than-one-pixel kind of effect.

Note also that if you change the array size from 128x128 -> 127x127, then the discrepancy is ~5 orders of magnitude smaller:

[cwwalter]

I've been looking for (and so far haven't found) a thread on the Rubin slack where I described a similar problem. I was mixing numpy and scipy forward and inverse FFTs and there were subtle differences in the location of the center of the coordinate system which were causing similar issues. I think this also supports Rachel and Josh's suggestions above.

[cwwalter]

@jecampagne Since you can look at it take a peek at:

https://lsstc.slack.com/archives/C2JPL8U78/p1569271022132800

From Fred M.

"Gotcha. One last tip for when you get back. The inconsistency that I mentioned is due to where the “center” is chosen for an array with an even number of pixels. I don’t remember the details but it basically comes down to the choice of choosing the center-left pixel or the center-right pixel as the center for a dimension that has an even number of pixels. Some of the numpy/scipy internal algorithms with FFTs use the center-left while others use the center-right, so their slicing algorithms are tricky to mimick."

Try using an odd size (say 129 x 129) and see if the problem disappears.

[jecampagne]

Thanks @rmandelb and @cwwalter. I have a version w/o pixel response and the effect 2 order of magnitude smaller than the difference between Galsim & scipy fftconvolve. So, I guess the center location is a good track to follow, I
will let you know the result of the galsim.hsm.FindAdaptiveMom asap.

[jecampagne]

I think @cwwalter you point to the right source of pb : below see the difference of Galsim vs scipy_ fftconvolve usign either a 128x128 image or a 129x129 image.

So, the source is identified +1, but I guess this numerical artifact can blow up in a another use case.
What do you recommend?

[jecampagne]

BTW, I have used the pmelchior scarlet fft.py set of routines mentionned in the #software-dev and I get the same dipole difference wrt to GalSim convolve (with a different sign)

[rmjarvis]

This is clearly just a convention difference in where to put the center of the result. IMO, GalSIm does it right, but as long as you know what to expect in terms of centering for the others, then they aren't really wrong either. I think this can be closed.

[cwwalter]

As Mike says, this is a convention.

I found this difference when trying to mix the FFT/Convolve routines between numpy and scipy. Galsim wasn't even involved. So, the main lesson I took away from that is that you need to use a consistent set of those routines from the same package. Fred and Peter wrote their own routines for Scarlet partly for these reasons (see that discussion I linked above).

If you know if is just the centering convention, and you still want to mix things, you can use the trick of picking odd sizes to remove the ambiguity. [I think trying to adjust things back and forth to make it work anyway is too dangerous, you probably won't be able to do it correctly consistently and is a recipe for mistakes]

[jecampagne]

Ok. close the thread. Sorry to have made troubles. mea culpa.

[cwwalter]

This was certainly confusing for me too when I first heard came across it and I also had to ask for help. So, I'm glad you were able to understand your use case too.