PyDL:docs

- Index
- Modules

### Navigation

- pydl v1.0.0 »
- Module code »
- pydl »
- pydl.photoop »

# 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

    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