Commit a0076109 authored by karpov-sv's avatar karpov-sv
Browse files

PSF documentation added

parent 1e61322d
=============================================
STDPipe - Simple Transient Detection Pipeline
---------------------------------------------
=============================================
*AKA: random codes noone else will ever use*
......@@ -26,10 +27,11 @@ User Guide
.. toctree::
:maxdepth: 3
About <self>
installation
usage
todo
contributing
todo
API documentation <api/modules>
Contributing
......
Point Spread Function (PSF) models
==================================
*STDPipe* includes basic support for point spread function (PSF) construction and analysis through the interface to `PSFEx <https://github.com/astromatic/psfex>`_ code that is able to build supersampled PSF models from the object lists created using `SExtractor <https://github.com/astromatic/sextractor>`_.
Please consider checking `its documentation <https://psfex.readthedocs.io>`__ to better understand the concepts of PSFEx operation and its possible configuration options.
The :func:`stdpipe.psf.run_psfex` routine transparently calls the `SExtractor`_ on the image and then uses its output (it needs a special set of measured object parameters and postage stamps for them) to run `PSFEx`_. The results are then parsed and returned as a structure representing the PSF model file as described `here <https://psfex.readthedocs.io/en/latest/Appendices.html#psf-file-format-description>`__. This structure may later be used to visualize PSF shape at arbitrary position inside the image, as well as to generate artificial stars and inject them into the image.
.. autofunction:: stdpipe.psf.run_psfex
:noindex:
.. attention::
While `PSFEx`_ may build PSF models depending on a variety of object parameters as described `here <https://psfex.readthedocs.io/en/latest/Working.html#managing-psf-variations>`__, *STDPipe* currently allows reading and interpreting only the ones depending on the position on the image (e.g. representing optical distortions). On the other hand, the original PSFEx model file may be stored to the file using `psffile` parameter of :func:`stdpipe.psf.run_psfex`, and then supplied directly to `SExtractor`_ through `psf` parameter of :func:`stdpipe.photometry.get_objects_sextractor` to be used for SExtractor PSF photometry.
Using PSF models
----------------
The PSF model built by `PSFEx`_ may also be loaded from PSFEx model file. This model may then be used to construct PSF stamps both in supersampled PSF model pixels, or original image pixels.
.. autofunction:: stdpipe.psf.load_psf
:noindex:
.. autofunction:: stdpipe.psf.get_supersampled_psf_stamp
:noindex:
.. autofunction:: stdpipe.psf.get_psf_stamp
:noindex:
Placing artificial stars into the image
---------------------------------------
*STDPipe* also contains a couple of convenience functions to directly place a PSF model at a given image position - that would correspond to artificial point source injection into the image. Higher-level one, :func:`stdpipe.pipeline.place_random_stars`, inserts a number of randomly generated sources with realistically (log-uniform in fluxes) distributed brightness at random positions, as well as returns the catalogue of injected stars that may be used e.g. in deriving the object detection efficiency.
.. code-block:: python
# Derive PSF model assuming that it does not change over the image
psf_model = psf.run_psfex(image, mask=mask, order=0, verbose=True)
# ..and now perform the detection efficiency analysis
sims = [] # will hold simulated stars
# We will repeatedly inject the stars (20 at a time), detect objects, and
# match them in order to see whether simulated stars are detectable
for _ in tqdm(range(100)):
# We will operate on a copy of original image
image1 = image.copy()
# We will use this value as a saturation level
saturation = 50000
# Simulate 20 random stars
sim = pipeline.place_random_stars(image1, psf_model, nstars=20, minflux=10, maxflux=1e6,
wcs=wcs, gain=gain, saturation=saturation)
# Initial metadata for injected stars
sim['mag_calib'] = sim['mag'] + zero_fn(sim['x'], sim['y']) # Apply zero point from photometric calibration
sim['detected'] = False
sim['mag_measured'] = np.nan
sim['magerr_measured'] = np.nan
sim['flags_measured'] = np.nan
# Mask corresponding to the saturation level
mask1 = image1 >= saturation
obj1 = photometry.get_objects_sextractor(image1, mask=mask|mask1, r0=1, aper=5.0,
wcs=wcs, gain=gain, minarea=3, sn=5)
# Apply zero point from photometric calibration
obj1['mag_calib'] = obj1['mag'] + zero_fn(obj1['x'], obj1['y'])
# Positional match within FWHM/2 radius
oidx,sidx,dist = astrometry.spherical_match(obj1['ra'], obj1['dec'], sim['ra'], sim['dec'],
pixscale*np.median(obj1['fwhm'])/2)
# Mark matched stars
sim['detected'][sidx] = True
# Also store measured magnitude, its error and flags
sim['mag_measured'][sidx] = obj1['mag_calib'][oidx]
sim['magerr_measured'][sidx] = obj1['magerr'][oidx]
sim['flags_measured'][sidx] = obj1['flags'][oidx]
sims.append(sim)
# Now we may stack the data from all runs into a single table
from astropy.table import vstack
sims = vstack(sims)
# ..and plot it as a histograms of magnitudes
h0,b0,_ = plt.hist(sims['mag_calib'], range=[12,22], bins=50, alpha=0.3, label='All simulated stars');
h1,b1,_ = plt.hist(sims['mag_calib'][sims['detected']], range=[12,22], bins=50, alpha=0.3, label='Detected');
h2,b2,_ = plt.hist(sims['mag_calib'][idx], range=[12,22], bins=50, alpha=0.3, label='Detected and unflagged');
plt.legend()
plt.xlabel('Injected magnitude')
plt.show()
# ..or as a detection efficiency
plt.plot(0.5*(b0[1:]+b0[:-1]), h1/h0, 'o-', label='Detected')
plt.plot(0.5*(b0[1:]+b0[:-1]), h2/h0, 'o-', label='Detected and unflagged')
plt.legend()
plt.xlabel('Injected magnitude')
plt.title('Fraction of detected artificial stars')
The complete example of detection efficiency analysis may be seen in the `corresponding notebook <https://github.com/karpov-sv/stdpipe/blob/master/notebooks/simulated_stars.ipynb>`_.
.. attention::
Note that the flux that you specify for the artificial star is a total flux - it will not necessarily correspond to the one measured inside an aperture (especially if the aperture size is quite small). The difference (*aperture correction*) is due to the fraction of the total flux that falls outside of your aperture, and it has to be somehow estimated and corrected if you wish to compare the generated fluxes of injected stars with measured values.
.. autofunction:: stdpipe.psf.place_psf_stamp
:noindex:
.. autofunction:: stdpipe.pipeline.place_random_stars
:noindex:
PSF photometry in SExtractor
----------------------------
Derived PSF models may also be used for performing PSF photometry in SExtractor.
As currently *STDPipe* lacks the routine for saving PSF models back to FITS files, in order to use PSF photometry you will need to pass `psffile` argument to :func:`stdpipe.psf.run_psfex` in order to directly store the produced PSF model into the file. Then, this file may be used in SExtractor by passing its path as a `psf` argument into :func:`stdpipe.photometry.get_objects_sextractor`. It will then produce the list of detected objects having a set of additional measured parameters including `flux_psf`, `fluxerr_psf`, `mag_psf`, `magerr_psf` - they correspond to the flux derived from fitting the objects with the PSF model.
.. code-block:: python
# We probably do not have enough stars to study PSF spatial variance, so we use order=0 here
psf_model = psf.run_psfex(image, mask=mask, order=0, gain=gain, psffile='/tmp/psf.psf', verbose=True)
obj_psf = photometry.get_objects_sextractor(image, mask=mask, aper=3.0, edge=10, gain=gain,
wcs=wcs, psf='/tmp/psf.psf', verbose=False)
# Now we may calibrate the photometry using PSF fluxes
m_psf = pipeline.calibrate_photometry(obj_psf, cat, sr=1/3600,
obj_col_mag='mag_psf', obj_col_mag_err='magerr_psf',
cat_col_mag='rmag', cat_col_mag1='gmag', cat_col_mag2='rmag',
order=0, verbose=True)
......@@ -48,6 +48,7 @@ Data processing
subtraction
transients
cutouts
psf
Common principles
-----------------
......
This diff is collapsed.
......@@ -346,14 +346,28 @@ def calibrate_photometry(obj, cat, sr=None, pixscale=None, order=0, bg_order=Non
return m
def make_random_stars(width=None, height=None, shape=None, nstars=100, minflux=1, maxflux=100000, gain=1, edge=0, wcs=None, verbose=False):
"""
Generate a table of random stars.
def make_random_stars(width=None, height=None, shape=None, nstars=100, minflux=1, maxflux=100000, edge=0, wcs=None, verbose=False):
"""Generate a table of random stars.
Coordinates are distributed uniformly with :code:`edge <= x < width-edge` and :code:`edge <= y < height-edge`.
Coordinates are distributed uniformly.
Fluxes are log-uniform between user-provided min and max values.
Returns: the catalogue of generated stars, with x, y and flux fields set.
Returns the catalogue of generated stars, with at least `x`, `y` and `flux` fields set.
If `wcs` is set, the returned catalogue will also contain `ra` and `dec` fields with
sky coordinates of the stars.
:param width: Width of the image containing generated stars
:param height: Height of the image containing generated stars
:param shape: Shape of the image containing generated stars, to be used instead of `width` and `height` if set
:param nstars: Number of artificial stars to inject into the image
:param minflux: Minimal flux of arttificial stars, in ADU units
:param maxflux: Maximal flux of arttificial stars, in ADU units
:param edge: Minimal distance to image edges for artificial stars to be placed. Optional
:param wcs: WCS as :class:`astropy.wcs.WCS` to be used to derive sky coordinates of injected stars. Optional
:param verbose: Whether to show verbose messages during the run of the function or not. May be either boolean, or a `print`-like function.
:returns: The catalogue of injected stars, containing the fluxes, image coordinates and (if `wcs` is set) sky coordinates of injected stars.
"""
# Simple wrapper around print for logging in verbose mode only
......@@ -378,20 +392,41 @@ def make_random_stars(width=None, height=None, shape=None, nstars=100, minflux=1
return cat
def place_random_stars(image, psf_model, nstars=100, minflux=1, maxflux=100000, gain=1, saturation=65535, edge=0, wcs=None, verbose=False):
"""
Randomly place artificial stars into the image.
Coordinates are distributed uniformly.
Fluxes are log-uniform between user-provided min and max values.
def place_random_stars(image, psf_model, nstars=100, minflux=1, maxflux=100000, gain=None, saturation=65535, edge=0, wcs=None, verbose=False):
"""Randomly place artificial stars into the image.
The stars will be placed on top of the existing content of the image, and the Poissonian
noise will be applied to these stars according to the specified `gain` value. Also, the saturation
level will be applied to the resulting image according to the `saturation` value.
Coordinates of the injected stars are distributed uniformly.
Fluxes are log-uniform between user-provided `minflux` and `maxflux` values in ADU units.
Returns the catalogue of generated stars, with at least `x`, `y` and `flux` fields set,
as returned by :func:`stdpipe.pipeline.make_random_stars`
If `wcs` is set, the returned catalogue will also contain `ra` and `dec` fields with
sky coordinates of the stars
:param image: Image where artificial stars will be injected
:param psf_model: PSF model structure as returned by :func:`stdpipe.psf.run_psfex`
:param nstars: Number of artificial stars to inject into the image
:param minflux: Minimal flux of arttificial stars, in ADU units
:param maxflux: Maximal flux of arttificial stars, in ADU units
:param gain: Image gain value. If set, will be used to apply Poissonian noise to the source
:param saturation: Saturation level in ADU units to be applied to the image
:param edge: Minimal distance to image edges for artificial stars to be placed
:param wcs: WCS as :class:`astropy.wcs.WCS` to be used to derive sky coordinates of injected stars
:param verbose: Whether to show verbose messages during the run of the function or not. May be either boolean, or a `print`-like function
:returns: The catalogue of injected stars, containing the fluxes, image coordinates and (if `wcs` is set) sky coordinates of injected stars.
Returns: the catalogue of generated stars, with x, y and flux fields set.
"""
# Simple wrapper around print for logging in verbose mode only
log = (verbose if callable(verbose) else print) if verbose else lambda *args,**kwargs: None
cat = make_random_stars(shape=image.shape, nstars=nstars, minflux=minflux, maxflux=maxflux,
gain=gain, edge=edge, wcs=wcs, verbose=verbose)
edge=edge, wcs=wcs, verbose=verbose)
for _ in cat:
psf.place_psf_stamp(image, psf_model, _['x'], _['y'], flux=_['flux'], gain=gain)
......
"""
Module for working with point-spread function (PSF) models
"""
from __future__ import absolute_import, division, print_function, unicode_literals
import os, shutil, tempfile, shlex
......@@ -9,12 +13,47 @@ from astropy.table import Table
from . import photometry
from . import utils
def run_psfex(image, mask=None, thresh=2.0, aper=3.0, r0=0.5, gain=1, minarea=5, vignet_size=None, order=0, sex_opts={}, checkimages=[], extra={}, psffile=None, _workdir=None, _tmpdir=None, _exe=None, verbose=False):
"""
Wrapper around PSFEx
def run_psfex(image, mask=None, thresh=2.0, aper=3.0, r0=0.5, gain=1, minarea=5, vignet_size=None, order=0, sex_extra={}, checkimages=[], extra={}, psffile=None, _workdir=None, _tmpdir=None, _exe=None, _sex_exe=None, verbose=False):
"""Wrapper around PSFEx to help extracting PSF models from images.
For the details of PSFEx operation we suggest to consult its documentation at https://psfex.readthedocs.io
:param image: Input image as a NumPy array
:param mask: Image mask as a boolean array (True values will be masked), optional
:param thresh: Detection threshold in sigmas above local background, for running initial SExtractor object detection
:param aper: Circular aperture radius in pixels, to be used for initial SExtractor object detection
:param r0: Smoothing kernel size (sigma) to be used for improving object detection in initial SExtractor call
:param gain: Image gain
:param minarea: Minimal number of pixels in the object to be considered a detection (`DETECT_MINAREA` parameter of SExtractor)
:param vignet_size: The size of *postage stamps* to be used for PSF model creation
:param order: Spatial order of PSF model variance
:param sex_extra: Dictionary of additional options to be passed to SExtractor for initial object detection (`extra` parameter of :func:`stdpipe.photometry.get_objects_sextractor`). Optional
:param checkimages: List of PSFEx checkimages to return along with PSF model. Optional.
:param extra: Dictionary of extra configuration parameters to be passed to PSFEx call, with keys as parameter names. See :code:`psfex -dd` for the full list.
:param psffile: If specified, PSF model file will also be stored under this file name, so that it may e.g. be re-used by SExtractor later. Optional
:param _workdir: If specified, all temporary files will be created in this directory, and will be kept intact after running SExtractor and PSFEx. May be used for debugging exact inputs and outputs of the executable. Optional
:param _tmpdir: If specified, all temporary files will be created in a dedicated directory (that will be deleted after running the executable) inside this path.
:param _exe: Full path to PSFEx executable. If not provided, the code tries to locate it automatically in your :envvar:`PATH`.
:param _sex_exe: Full path to SExtractor executable. If not provided, the code tries to locate it automatically in your :envvar:`PATH`.
:param verbose: Whether to show verbose messages during the run of the function or not. May be either boolean, or a `print`-like function.
:returns: PSF structure corresponding to the built PSFEx model.
The structure has at least the following fields:
- `width`, `height` - dimesnions of supersampled PSF stamp
- `fwhm` - mean full width at half maximum (FWHM) of the images used for building the PSF model
- `sampling` - conversion factor between PSF stamp (supersampled) pixel size, and original image one (less than unity when supersampled resolution is finer than original image one)
- `ncoeffs` - number of coefficients pixel polynomials have
- `degree` - polynomial degree of a spatial variance of PSF model
- `data` - the data containing per-pixel polynomial coefficients for PSF model
- `header` - original FITS header of PSF model file, if :code:`get_header=True` parameter was set
This structure corresponds to the contents of original PSFEx generated output file that
is documented at https://psfex.readthedocs.io/en/latest/Appendices.html#psf-file-format-description
"""
# Simple wrapper around print for logging in verbose mode only
# Simple wrapper around print for logging in verbose mode only
log = (verbose if callable(verbose) else print) if verbose else lambda *args,**kwargs: None
# Find the binary
......@@ -45,7 +84,7 @@ def run_psfex(image, mask=None, thresh=2.0, aper=3.0, r0=0.5, gain=1, minarea=5,
log('Extracting PSF using vignette size %d x %d pixels' % (vignet_size, vignet_size))
# Run SExtractor on input image in current workdir so that the LDAC catalogue will be in out.cat there
obj = photometry.get_objects_sextractor(image, mask=mask, thresh=thresh, aper=aper, r0=r0, gain=gain, minarea=minarea, _workdir=workdir, _tmpdir=_tmpdir, verbose=verbose, extra_params=['SNR_WIN', 'ELONGATION', 'VIGNET(%d,%d)' % (vignet_size,vignet_size)], extra=sex_opts)
obj = photometry.get_objects_sextractor(image, mask=mask, thresh=thresh, aper=aper, r0=r0, gain=gain, minarea=minarea, _workdir=workdir, _tmpdir=_tmpdir, _exe=_sex_exe, verbose=verbose, extra_params=['SNR_WIN', 'ELONGATION', 'VIGNET(%d,%d)' % (vignet_size,vignet_size)], extra=sex_extra)
catname = os.path.join(workdir, 'out.cat')
psfname = os.path.join(workdir, 'out.psf')
......@@ -108,11 +147,18 @@ def run_psfex(image, mask=None, thresh=2.0, aper=3.0, r0=0.5, gain=1, minarea=5,
return result
def load_psf(filename, get_header=False, verbose=False):
"""
Load PSF model from PSFEx file
"""Load PSF model from PSFEx file
The structure may be useful for inspection of PSF model with :func:`stdpipe.psf.get_supersampled_psf_stamp` and :func:`stdpipe.psf.get_psf_stamp`, as well as for injection of PSF instances (fake objects) into the image with :func:`stdpipe.psf.place_psf_stamp`.
:param filename: Name of a file containing PSF model built by PSFEx
:param get_header: Whether to return the original FITS header of PSF model file or not. If set, the header will be stored in `header` field of the returned structire
:param verbose: Whether to show verbose messages during the run of the function or not. May be either boolean, or a `print`-like function.
:returns: PSF structure in the same format as returned from :func:`stdpipe.psf.run_psfex`.
"""
# Simple wrapper around print for logging in verbose mode only
# Simple wrapper around print for logging in verbose mode only
log = (verbose if callable(verbose) else print) if verbose else lambda *args,**kwargs: None
log('Loading PSF model from %s' % filename)
......@@ -175,8 +221,18 @@ def bilinear_interpolate(im, x, y):
return wa*Ia + wb*Ib + wc*Ic + wd*Id
def get_supersampled_psf_stamp(psf, x=0, y=0, normalize=True):
"""
Returns supersampled PSF model for a given image position
"""Returns supersampled PSF model for a given position inside the image.
The returned stamp corresponds to PSF model evaluated at a given position inside the image,
with its center always in the center of central stamp pixel.
Every *supersampled* pixel of the stamp corresponds to :code:`psf['sampling']` pixels of the original image.
:param psf: Input PSF structure as returned by :func:`stdpipe.psf.run_psfex` or :func:`stdpipe.psf.load_psf`
:param x: `x` coordinate of the position inside the original image to evaluate the PSF model
:param y: `y` coordinate of the position inside the original image to evaluate the PSF model
:param normalize: Whether to normalize the stamp to have flux exactly equal to unity or not
:returns: stamp of the PSF model evaluated at the given position inside the image
"""
dx = 1.0*(x - psf['x0'])/psf['sx']
......@@ -196,18 +252,34 @@ def get_supersampled_psf_stamp(psf, x=0, y=0, normalize=True):
return stamp
def get_psf_stamp(psf, x=0, y=0, dx=None, dy=None, normalize=True):
"""
Returns PSF stamp in image space with sub-pixel shift applied.
Stamp is odd-sized, with center at
"""Returns PSF stamp in original image pixel space with sub-pixel shift applied.
The PSF model is evaluated at the requested position inside the original image,
and then downscaled from supersampled pixels of the PSF model to original image pixels,
and then adjusted to accommodate for requested :code:`(dx, dy)` sub-pixel shift.
Stamp is odd-sized, with PSF center at::
x0 = floor(width/2) + dx
y0 = floor(height/2) + dy
If dx or dy is None, they are computed as
If :code:`dx=None` or :code:`dy=None`, they are computed directly from the
floating point parts of the position `x` and `y` arguments::
dx = x - round(x)
dy = y - round(y)
The stamp should directly represent stellar shape at a given position (including sub-pixel
center shift) inside the image.
:param psf: Input PSF structure as returned by :func:`stdpipe.psf.run_psfex` or :func:`stdpipe.psf.load_psf`
:param x: `x` coordinate of the position inside the original image to evaluate the PSF model
:param y: `y` coordinate of the position inside the original image to evaluate the PSF model
:param dx: Sub-pixel adjustment of PSF position in image space, `x` direction
:param dy: Sub-pixel adjustment of PSF position in image space, `y` direction
:param normalize: Whether to normalize the stamp to have flux exactly equal to unity or not
:returns: Stamp of the PSF model evaluated at the given position inside the image, in original image pixels.
"""
if dx is None:
......@@ -247,10 +319,24 @@ def get_psf_stamp(psf, x=0, y=0, dx=None, dy=None, normalize=True):
return stamp
def place_psf_stamp(image, psf, x0, y0, flux=1, gain=None):
"""
Places PSF stamp, scaled to a given flux, at a given position inside the image.
The stamp values are added to current content of the image.
If gain value is set, the Poissonian noise is applied to the stamp.
"""Places PSF stamp, scaled to a given flux, at a given position inside the image.
PSF stamp is evaluated at a given position, then adjusted to accommodate for
required sub-pixel shift, and finally scaled to requested flux value. Thus,
the routine corresponds to injection of an artificial point source into the image.
The stamp values are added on top of current content of the image.
If `gain` value is set, the Poissonian noise is applied to the stamp.
The image is modified in-place.
:param image: The image where artifi
:param psf: Input PSF structure as returned by :func:`stdpipe.psf.run_psfex` or :func:`stdpipe.psf.load_psf`
:param x0: `x` coordinate of the position to inject the source
:param y0: `y` coordinate of the position to inject the source
:param flux: The source flux in ADU units
:param gain: Image gain value. If set, used to apply Poissonian noise to the source.
"""
stamp = get_psf_stamp(psf, x0, y0, normalize=True)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment