Tutorial - Part #2 - Spatially variant PSF

This package hosts a simple implementation of the proposed technique by Lauer 2002.

The technique performs a linear decomposition of the PSF into several basis elements of small size (much smaller than the image) normalized to unit sum. Each of these autopsfs come with a coefficient, which is not normalized, that have the same size of the image.

To perform the decomposition uses the Karhunen-Loeve transform, in order to calculate the basis as a linear expansion of the psf observations –which are cutout stamps of the stars in the image.

import numpy as np
import matplotlib.pyplot as plt

from properimage import single_image as s

%matplotlib inline

import os
import shlex
import subprocess
import sys

import numpy as np
import matplotlib.pyplot as plt

from properimage import propercoadd as pc
from properimage import single_image as si
from properimage import utils
from properimage import plot
from properimage import simtools

We simulate an image with several star shapes across the field

frames = []
for theta in [0, 45, 105, 150]:
    N = 512  # side
    X_FWHM = 5 + 5.5*theta/180
    Y_FWHM = 5
    t_exp = 5
    max_fw = max(X_FWHM, Y_FWHM)
    #test_dir = os.path.abspath('./test/test_images/psf_basis_kl_gs')

    x = np.random.randint(low=6*max_fw, high=N-6*max_fw, size=80)
    y = np.random.randint(low=6*max_fw, high=N-6*max_fw, size=80)
    xy = [(x[i], y[i]) for i in range(80)]

    SN =  30. # SN para poder medir psf
    weights = list(np.linspace(10, 1000., len(xy)))
    m = simtools.delta_point(N, center=False, xy=xy, weights=weights)
    im = simtools.image(m, N, t_exp, X_FWHM, Y_FWHM=Y_FWHM, theta=theta,
                        SN=SN, bkg_pdf='gaussian')
    frames.append(im+100.)

mean = 0, std = 1.1570077460537824, b = 34.71023238161347, SN = 30.0
mean = 0, std = 0.8884857113082378, b = 26.654571339247134, SN = 30.0
mean = 0, std = 0.6862562794039239, b = 20.587688382117715, SN = 30.0
mean = 0, std = 0.6398376452574939, b = 19.195129357724817, SN = 30.0

frame = np.zeros((1024, 1024))
for j in range(2):
    for i in range(2):
        frame[i*512:(i+1)*512, j*512:(j+1)*512] = frames[i+2*j]

plt.figure(figsize=(12, 12))
plt.imshow(np.log(frame), interpolation='none')
plt.colorbar(shrink=0.7)

<matplotlib.colorbar.Colorbar at 0x7fe782f613a0>

We perform the psf extraction. This is possible to do with a context manager, in order to autoclean the disk files where the star stamps are stored.

with si.SingleImage(frame, smooth_psf=False) as sim:
    a_fields, psf_basis = sim.get_variable_psf(inf_loss=0.15)
    x, y = sim.get_afield_domain()

updating stamp shape to (25,25)
(109, 109) (625, 109)
cleaning...

print(len(psf_basis))

6

If we plot the coefficients, they look as smooth fields of the same size of the image.

Whenever they grow they show where the corresponding autopsf element is more important than others.

plot.plot_afields(a_fields=a_fields, x=x, y=y, nbook=True)

The autopsf elements are small patches, like the ones below.

The principal component psf is the first element of this list of arrays.

plot.plot_psfbasis(psf_basis=psf_basis, nbook=True)