Source code for pydl.photoop.image

# Licensed under a 3-clause BSD style license - see LICENSE.rst
# -*- coding: utf-8 -*-
"""This module corresponds to the image directory in photoop.
"""
import numpy as np


[docs]def sdss_psf_recon(psfield, xpos, ypos, normalize=None, trimdim=None): """Reconstruct the PSF at position (`xpos`, `ypos`) in an SDSS field. These can be read from the `psField files`_. .. _`psField files`: https://data.sdss.org/datamodel/files/PHOTO_REDUX/RERUN/RUN/objcs/CAMCOL/psField.html Parameters ---------- psfield : :class:`~astropy.io.fits.fitsrec.FITS_rec` A PSF KL-decomposition structure read from a psField file. This is from HDU's 1 through 5 in a psField file, corresponding to the five filters *u*, *g*, *r*, *i*, *z*. xpos : :class:`int` Column position (0-indexed, not 0.5-indexed as PHOTO outputs). ypos : :class:`int` Row position (0-indexed, not 0.5-indexed as PHOTO outputs). normalize : :class:`int` or :class:`float`, optional If set, normalize the integral of the image to this value. trimdim : :class:`tuple`, optional Trimmed dimensions; for example, set to (25, 25) to trim the output PSF image to those dimensions. These dimensions must be odd-valued. Returns ------- :class:`numpy.ndarray` The 2D reconstructed PSF image, typically dimensioned (51, 51). The center of the PSF is always the central pixel; this function will not apply any sub-pixel shifts. Notes ----- The SDSS photo PSF is described as a set of eigen-templates, where the mix of these eigen-templates is a simple polynomial function with (x,y) position on the CCD. Typically, there are 4 such 51x51 pixel templates. The polynomial functions are typically quadratic in each dimension, with no cross-terms. The formula is the following, where :math:`i` is the index of row polynomial order, :math:`j` is the index of column polynomial order, and :math:`k` is the template index: .. math:: a_k &= \\sum_i \\sum_j (0.001 \\times \\mathrm{ROWC})^i \\times (0.001 \\times \\mathrm{COLC})^j \\times [C_k]_{ij} \\\\ \\mathrm{psfimage} &= \\sum_k a_k \\mathrm{RROWS}_k The polynomial terms need not be of the same order for each template. Examples -------- >>> import numpy as np >>> from astropy.io import fits >>> from astropy.utils.data import get_readable_fileobj >>> from pydl.photoop.image import sdss_psf_recon >>> psfile = ('https://data.sdss.org/sas/dr14/eboss/photo/redux/301/' + ... '3366/objcs/3/psField-003366-3-0110.fit') >>> with get_readable_fileobj(psfile, encoding='binary') as psField: # doctest: +REMOTE_DATA ... with fits.open(psField) as hdulist: ... psf = sdss_psf_recon(hdulist[3].data, 500, 500) """ # # Hard-wired scale factor for ypos/xpos coefficients. # rc_scale = 0.001 # # Assume that the dimensions of each eigen-template are the same. # rnrow = psfield['RNROW'][0] rncol = psfield['RNCOL'][0] npix = rnrow * rncol # # These are the polynomial coefficients as a function of x,y. # Only compute these coefficients for the maximum polynomial order in use. # In general, this order can be different for each eigen-template. # nr_max = psfield['nrow_b'].max() nc_max = psfield['ncol_b'].max() # nb = nrow_b * ncol_b coeffs = np.outer(((ypos+0.5) * rc_scale)**np.arange(nr_max), ((xpos+0.5) * rc_scale)**np.arange(nc_max)) # # Reconstruct the image by summing each eigen-template. # neigen = psfield.shape[0] psfimage = np.zeros(psfield['RROWS'][0].shape, dtype=psfield['RROWS'][0].dtype) for i in range(neigen): psfimage += (psfield['RROWS'][i] * (psfield[i]['c'][0:psfield[0]['nrow_b'], 0:psfield[0]['ncol_b']] * coeffs[0:psfield[0]['nrow_b'], 0:psfield[0]['ncol_b']].T).sum()) # # We have reconstructed the PSF as a vector using all the # pixels in psfield.RROWS. So, at the end we trim this vector to only # those pixels used, and reform it into a 2-dimensional image. # if normalize is not None: psfimage /= psfimage.sum() psfimage *= normalize psfimage = psfimage[0:npix].reshape(rncol, rnrow) if trimdim is not None: psfimage = psfimage[(rncol - trimdim[0])//2:(rncol + trimdim[0])//2, (rnrow - trimdim[1])//2:(rnrow + trimdim[1])//2] return psfimage