# 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) 