Add function to produce metacal images

[EiffL]

I've added a proof of concept notebook here to produce metacal images, i.e. images that have an artificial shear added to them. This notebook is doing everything in Fourier space (deconvolution, shear, reconvolution) before returning an image in pixel space. And currently assumes that all input quantities are provided as Fourier arrays.

So, we want a generic function, that would look like something like this:

from autometacal.metacal import generate_mcal_images

# [...]

metacal_image = generate_mcal_images(galaxy_stamps, noise_stamps, PSF, g)

where g is the added shear tensor. The noise is a noise realisation following the statistics of the noise in the image, that will be sheared with negative shear to cancel out the effect of shearing the noise.

[EiffL]

Here is an example of measuring moments with ngmix, we can probably use something like that to measure shear on the metacal images produced in this issue. https://github.com/esheldon/ngmix/blob/master/ngmix/moments.py

[aguinot]

Here are the step to handle the metacal process in ngmix:

• Setup data:
  • https://github.com/esheldon/ngmix/blob/a0642277eadc86c80adb560d237c1cada82c12f7/ngmix/metacal.py#L629
  • https://github.com/esheldon/ngmix/blob/a0642277eadc86c80adb560d237c1cada82c12f7/ngmix/metacal.py#L755
• Once the data are setup, it is trivial to apply the metacal. It comes down to call .shear to apply the shear and then .convolve with the new PSF.

I would recommend to have a clear setup step which could even happen before going through the autodiff. Here is a suggestion for the processing (I think it is similar to what Francois has done):

• Create an interpolated image in Fourier space deconvolve from the original PSF.
• Create the new PSF as a galsim.GSObject or galsim.InterpolatedImage (depending which strategy is used for the new PSF: gaussian, gaussian mixture, symmetrization,...). And make the Fourier transform.
• Create a galsim.WCSBase which contain the local jacobian for the WCS.
• Make the autodiff where the argument would be the interpolated image, the new PSF, the WCS (I think that is all we need).
• The autodiff only does shearing, convolution and shape measurement. All those step can be done in Fourier.

[EiffL]

I'm going to work on this in branch u/EiffL/mcal_imags

[andrevitorelli]

I am currently working on a function to get the psf deconvolution kernel as a tf tensor. So, we get an image as a PSF model, convert them to a galsim object, apply the Deconvolve method, and convert to the output. I guess my main question is about image size. We keep the size of the image as the original psf image size? How does that relate to the galaxy image stamp size?

An example of the code to illustrate the idea:

def PSFdeconvolve(psf_img,
                  pixel_scale=0.2,
                  interp_factor=2,
                  padding_factor=1):
    """
    Returns a deconvolution kernel of a psf image.

    Args:
        psf_img: numpy array representing the psf model
        pixel_scale: the pixel scale of the image, in arcsec/pixel
        interp_factor: the interpolation factor for super-resolution
        padding_factor: a factor to add side pads to the image
    Returns:
        A complex tensorflow tensor that is a deconvolution kernel.

    """

    N = len(psf_img)

    psf_galsim=galsim.InterpolatedImage(galsim.Image(psf_img,scale=pixel_scale))

    ipsf=galsim.Deconvolve(psf_galsim)

    Nk = N*interp_factor*padding_factor
    bounds = galsim._BoundsI(-Nk//2, Nk//2-1, -Nk//2, Nk//2-1)

    imipsf = ipsf.drawKImage(bounds=bounds, scale=2.*np.pi/(N*padding_factor* pixel_scale), recenter=False)

    return tf.convert_to_tensor(imipsf.array,dtype=tf.complex64)

[EiffL]

Annnd with the merging of #16 we have a first prototype of this, I'm going to close this issue for now