# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
This module provides helper utilities for making PSF models.
"""
import contextlib
import re
import numpy as np
from astropy.modeling import CompoundModel, Fittable2DModel, Parameter
from astropy.modeling.models import Const2D, Identity, Shift
from astropy.nddata import NDData
from astropy.units import Quantity
from astropy.utils.decorators import deprecated
from scipy.integrate import dblquad, trapezoid
__all__ = ['PRFAdapter', 'grid_from_epsfs', 'make_psf_model']
[docs]
def make_psf_model(model, *, x_name=None, y_name=None, flux_name=None,
normalize=True, dx=50, dy=50, subsample=100,
use_dblquad=False):
"""
Make a PSF model that can be used with the PSF photometry classes
(`PSFPhotometry` or `IterativePSFPhotometry`) from an Astropy
fittable 2D model.
If the ``x_name``, ``y_name``, or ``flux_name`` keywords are input,
this function will map those ``model`` parameter names to ``x_0``,
``y_0``, or ``flux``, respectively.
If any of the ``x_name``, ``y_name``, or ``flux_name`` keywords
are `None`, then a new parameter will be added to the model
corresponding to the missing parameter. Any new position parameters
will be set to a default value of 0, and any new flux parameter will
be set to a default value of 1.
The output PSF model will have ``x_name``, ``y_name``, and
``flux_name`` attributes that contain the name of the corresponding
model parameter.
.. note::
This function is needed only in cases where the 2D PSF model
does not have ``x_0``, ``y_0``, and ``flux`` parameters.
It is *not* needed for any of the PSF models provided by
Photutils.
Parameters
----------
model : `~astropy.modeling.Fittable2DModel`
An Astropy fittable 2D model to use as a PSF.
x_name : `str` or `None`, optional
The name of the ``model`` parameter that corresponds to the x
center of the PSF. If `None`, the model will be assumed to be
centered at x=0, and a new model parameter called ``xpos_0``
will be added for the x position.
y_name : `str` or `None`, optional
The name of the ``model`` parameter that corresponds to the
y center of the PSF. If `None`, the model will be assumed
to be centered at y=0, and a new parameter called ``ypos_1``
will be added for the y position.
flux_name : `str` or `None`, optional
The name of the ``model`` parameter that corresponds to the
total flux of a source. If `None`, a new model parameter called
``flux_3`` will be added for model flux.
normalize : bool, optional
If `True`, the input ``model`` will be integrated and rescaled
so that its sum integrates to 1. This normalization occurs only
once for the input ``model``. If the total flux of ``model``
somehow depends on (x, y) position, then one will need to
correct the fitted model fluxes for this effect.
dx, dy : odd int, optional
The size of the integration grid in x and y for normalization.
Must be odd. These keywords are ignored if ``normalize`` is
`False` or ``use_dblquad`` is `True`.
subsample : int, optional
The subsampling factor for the integration grid along each axis
for normalization. Each pixel will be sampled ``subsample`` x
``subsample`` times. This keyword is ignored if ``normalize`` is
`False` or ``use_dblquad`` is `True`.
use_dblquad : bool, optional
If `True`, then use `scipy.integrate.dblquad` to integrate the
model for normalization. This is *much* slower than the default
integration of the evaluated model, but it is more accurate.
This keyword is ignored if ``normalize`` is `False`.
Returns
-------
result : `~astropy.modeling.CompoundModel`
A PSF model that can be used with the PSF photometry classes.
The returned model will always be an Astropy compound model.
Notes
-----
To normalize the model, by default it is discretized on a grid of
size ``dx`` x ``dy`` from the model center with a subsampling factor
of ``subsample``. The model is then integrated over the grid using
trapezoidal integration.
If the ``use_dblquad`` keyword is set to `True`, then the model is
integrated using `scipy.integrate.dblquad`. This is *much* slower
than the default integration of the evaluated model, but it is more
accurate. Also, note that the ``dblquad`` integration can sometimes
fail, e.g., return zero for a non-zero model. This can happen when
the model function is sharply localized relative to the size of the
integration interval.
Examples
--------
>>> from astropy.modeling.models import Gaussian2D
>>> from photutils.psf import make_psf_model
>>> model = Gaussian2D(x_stddev=2, y_stddev=2)
>>> psf_model = make_psf_model(model, x_name='x_mean', y_name='y_mean')
>>> print(psf_model.param_names) # doctest: +SKIP
('amplitude_2', 'x_mean_2', 'y_mean_2', 'x_stddev_2', 'y_stddev_2',
'theta_2', 'amplitude_3', 'amplitude_4')
"""
input_model = model.copy()
if x_name is None:
x_model = _InverseShift(0, name='x_position')
# "offset" is the _InverseShift parameter name;
# the x inverse shift model is always the first submodel
x_name = 'offset_0'
else:
if x_name not in input_model.param_names:
raise ValueError(f'{x_name!r} parameter name not found in the '
'input model.')
x_model = Identity(1)
x_name = _shift_model_param(input_model, x_name, shift=2)
if y_name is None:
y_model = _InverseShift(0, name='y_position')
# "offset" is the _InverseShift parameter name;
# the y inverse shift model is always the second submodel
y_name = 'offset_1'
else:
if y_name not in input_model.param_names:
raise ValueError(f'{y_name!r} parameter name not found in the '
'input model.')
y_model = Identity(1)
y_name = _shift_model_param(input_model, y_name, shift=2)
x_model.fittable = True
y_model.fittable = True
psf_model = (x_model & y_model) | input_model
if flux_name is None:
psf_model *= Const2D(1.0, name='flux')
# "amplitude" is the Const2D parameter name;
# the flux scaling is always the last component
flux_name = psf_model.param_names[-1]
else:
flux_name = _shift_model_param(input_model, flux_name, shift=2)
if normalize:
integral = _integrate_model(psf_model, x_name=x_name, y_name=y_name,
dx=dx, dy=dy, subsample=subsample,
use_dblquad=use_dblquad)
if integral == 0:
raise ValueError('Cannot normalize the model because the '
'integrated flux is zero.')
psf_model *= Const2D(1.0 / integral, name='normalization_scaling')
# fix all the output model parameters that are not x, y, or flux
for name in psf_model.param_names:
psf_model.fixed[name] = name not in (x_name, y_name, flux_name)
# final check that the x, y, and flux parameter names are in the
# output model
names = (x_name, y_name, flux_name)
for name in names:
if name not in psf_model.param_names:
raise ValueError(f'{name!r} parameter name not found in the '
'output model.')
# set the parameter names for the PSF photometry classes
psf_model.x_name = x_name
psf_model.y_name = y_name
psf_model.flux_name = flux_name
# set aliases
psf_model.x_0 = getattr(psf_model, x_name)
psf_model.y_0 = getattr(psf_model, y_name)
psf_model.flux = getattr(psf_model, flux_name)
return psf_model
class _InverseShift(Shift):
"""
A model that is the inverse of the normal
`astropy.modeling.functional_models.Shift` model.
"""
@staticmethod
def evaluate(x, offset):
return x - offset
@staticmethod
def fit_deriv(x, offset):
"""
One dimensional Shift model derivative with respect to
parameter.
"""
d_offset = -np.ones_like(x) + offset * 0.0
return [d_offset]
def _integrate_model(model, x_name=None, y_name=None, dx=50, dy=50,
subsample=100, use_dblquad=False):
"""
Integrate a model over a 2D grid.
By default, the model is discretized on a grid of size ``dx``
x ``dy`` from the model center with a subsampling factor of
``subsample``. The model is then integrated over the grid using
trapezoidal integration.
If the ``use_dblquad`` keyword is set to `True`, then the model is
integrated using `scipy.integrate.dblquad`. This is *much* slower
than the default integration of the evaluated model, but it is more
accurate. Also, note that the ``dblquad`` integration can sometimes
fail, e.g., return zero for a non-zero model. This can happen when
the model function is sharply localized relative to the size of the
integration interval.
Parameters
----------
model : `~astropy.modeling.Fittable2DModel`
The Astropy 2D model.
x_name : str or `None`, optional
The name of the ``model`` parameter that corresponds to the
x-axis center of the PSF. This parameter is required if
``use_dblquad`` is `False` and ignored if ``use_dblquad`` is
`True`.
y_name : str or `None`, optional
The name of the ``model`` parameter that corresponds to the
y-axis center of the PSF. This parameter is required if
``use_dblquad`` is `False` and ignored if ``use_dblquad`` is
`True`.
dx, dy : odd int, optional
The size of the integration grid in x and y. Must be odd.
These keywords are ignored if ``use_dblquad`` is `True`.
subsample : int, optional
The subsampling factor for the integration grid along each axis.
Each pixel will be sampled ``subsample`` x ``subsample`` times.
This keyword is ignored if ``use_dblquad`` is `True`.
use_dblquad : bool, optional
If `True`, then use `scipy.integrate.dblquad` to integrate the
model. This is *much* slower than the default integration of
the evaluated model, but it is more accurate.
Returns
-------
integral : float
The integral of the model over the 2D grid.
"""
if use_dblquad:
return dblquad(model, -np.inf, np.inf, -np.inf, np.inf)[0]
if dx <= 0 or dy <= 0:
raise ValueError('dx and dy must be > 0')
if subsample < 1:
raise ValueError('subsample must be >= 1')
xc = getattr(model, x_name)
yc = getattr(model, y_name)
if np.any(~np.isfinite((xc.value, yc.value))):
raise ValueError('model x and y positions must be finite')
hx = (dx - 1) / 2
hy = (dy - 1) / 2
nxpts = int(dx * subsample)
nypts = int(dy * subsample)
xvals = np.linspace(xc - hx, xc + hx, nxpts)
yvals = np.linspace(yc - hy, yc + hy, nypts)
# evaluate the model on the subsampled grid
data = model(xvals.reshape(-1, 1), yvals.reshape(1, -1))
if isinstance(data, Quantity):
data = data.value
# now integrate over the subsampled grid (first over x, then over y)
int_func = trapezoid
return int_func([int_func(row, xvals) for row in data], yvals)
def _shift_model_param(model, param_name, shift=2):
if isinstance(model, CompoundModel):
# for CompoundModel, add "shift" to the parameter suffix
out = re.search(r'(.*)_([\d]*)$', param_name)
new_name = out.groups()[0] + '_' + str(int(out.groups()[1]) + 2)
else:
# simply add the shift to the parameter name
new_name = param_name + '_' + str(shift)
return new_name
[docs]
def grid_from_epsfs(epsfs, grid_xypos=None, meta=None):
"""
Create a GriddedPSFModel from a list of ImagePSF models.
Given a list of `~photutils.psf.ImagePSF` models, this function will
return a `~photutils.psf.GriddedPSFModel`. The fiducial points for
each input ImagePSF can either be set on each individual model by
setting the 'x_0' and 'y_0' attributes, or provided as a list of
tuples (``grid_xypos``). If a ``grid_xypos`` list is provided, it
must match the length of input EPSFs. In either case, the fiducial
points must be on a grid.
Optionally, a ``meta`` dictionary may be provided for the output
GriddedPSFModel. If this dictionary contains the keys 'grid_xypos',
'oversampling', or 'fill_value', they will be overridden.
Note: If set on the input ImagePSF (x_0, y_0), then ``origin``
must be the same for each input EPSF. Additionally data units and
dimensions must be for each input EPSF, and values for ``flux`` and
``oversampling``, and ``fill_value`` must match as well.
Parameters
----------
epsfs : list of `photutils.psf.ImagePSF`
A list of ImagePSF models representing the individual PSFs.
grid_xypos : list, optional
A list of fiducial points (x_0, y_0) for each PSF. If not
provided, the x_0 and y_0 of each input EPSF will be considered
the fiducial point for that PSF. Default is None.
meta : dict, optional
Additional metadata for the GriddedPSFModel. Note that, if
they exist in the supplied ``meta``, any values under the keys
``grid_xypos`` , ``oversampling``, or ``fill_value`` will be
overridden. Default is None.
Returns
-------
GriddedPSFModel: `photutils.psf.GriddedPSFModel`
The gridded PSF model created from the input EPSFs.
"""
# prevent circular imports
from photutils.psf import GriddedPSFModel, ImagePSF
# optional, to store fiducial from input if `grid_xypos` is None
x_0s = []
y_0s = []
data_arrs = []
oversampling = None
fill_value = None
dat_unit = None
origin = None
flux = None
# make sure, if provided, that ``grid_xypos`` is the same length as
# ``epsfs``
if grid_xypos is not None and len(grid_xypos) != len(epsfs):
raise ValueError('``grid_xypos`` must be the same length as '
'``epsfs``.')
# loop over input once
for i, epsf in enumerate(epsfs):
# check input type
if not isinstance(epsf, ImagePSF):
raise TypeError('All input `epsfs` must be of type ImagePSF')
# get data array from EPSF
data_arrs.append(epsf.data)
if i == 0:
oversampling = epsf.oversampling
# same for fill value and flux, grid will have a single value
# so it should be the same for all input, and error if not.
fill_value = epsf.fill_value
# check that origins are the same
if grid_xypos is None:
origin = epsf.origin
flux = epsf.flux
# if there's a unit, those should also all be the same
with contextlib.suppress(AttributeError):
dat_unit = epsf.data.unit
else:
if np.any(epsf.oversampling != oversampling):
raise ValueError('All input ImagePSF models must have the '
'same value for ``oversampling``.')
if epsf.fill_value != fill_value:
raise ValueError('All input ImagePSF models must have the '
'same value for ``fill_value``.')
if epsf.data.ndim != data_arrs[0].ndim:
raise ValueError('All input ImagePSF models must have data '
'with the same dimensions.')
try:
unitt = epsf.data_unit
if unitt != dat_unit:
raise ValueError('All input data must have the same unit.')
except AttributeError as exc:
if dat_unit is not None:
raise ValueError('All input data must have the same '
'unit.') from exc
if epsf.flux != flux:
raise ValueError('All input ImagePSF models must have the '
'same value for ``flux``.')
if grid_xypos is None: # get gridxy_pos from x_0, y_0 if not provided
x_0s.append(epsf.x_0.value)
y_0s.append(epsf.y_0.value)
# also check that origin is the same, if using x_0s and y_0s
# from input
if np.all(epsf.origin != origin):
raise ValueError('If using ``x_0``, ``y_0`` as fiducial point,'
'``origin`` must match for each input EPSF.')
# if not supplied, use from x_0, y_0 of input EPSFs as fiducuals
# these are checked when GriddedPSFModel is created to make sure they
# are actually on a grid.
if grid_xypos is None:
grid_xypos = list(zip(x_0s, y_0s, strict=True))
data_cube = np.stack(data_arrs, axis=0)
if meta is None:
meta = {}
# add required keywords to meta
meta['grid_xypos'] = grid_xypos
meta['oversampling'] = oversampling
meta['fill_value'] = fill_value
data = NDData(data_cube, meta=meta)
return GriddedPSFModel(data, fill_value=fill_value)
[docs]
@deprecated('2.0.0', alternative='a FittableImageModel derived from the '
'discretize_model function in astropy.convolution')
class PRFAdapter(Fittable2DModel):
"""
A model that adapts a supplied PSF model to act as a PRF.
It integrates the PSF model over pixel "boxes". A critical built-in
assumption is that the PSF model scale and location parameters are
in *pixel* units.
Parameters
----------
psfmodel : a 2D model
The model to assume as representative of the PSF.
renormalize_psf : bool, optional
If True, the model will be integrated from -inf to inf and
rescaled so that the total integrates to 1. Note that this
renormalization only occurs *once*, so if the total flux of
``psfmodel`` depends on position, this will *not* be correct.
flux : float, optional
The total flux of the star.
x_0 : float, optional
The x position of the star.
y_0 : float, optional
The y position of the star.
xname : str or None, optional
The name of the ``psfmodel`` parameter that corresponds to the
x-axis center of the PSF. If None, the model will be assumed to
be centered at x=0.
yname : str or None, optional
The name of the ``psfmodel`` parameter that corresponds to the
y-axis center of the PSF. If None, the model will be assumed to
be centered at y=0.
fluxname : str or None, optional
The name of the ``psfmodel`` parameter that corresponds to
the total flux of the star. If None, a scaling factor will
be applied by the ``PRFAdapter`` instead of modifying the
``psfmodel``.
**kwargs : dict, optional
Additional optional keyword arguments to be passed to the
`astropy.modeling.Model` parent class.
Notes
-----
This current implementation of this class (using numerical
integration for each pixel) is extremely slow, and only suited for
experimentation over relatively few small regions. It should be
used only when absolutely necessary. If a model class of this type
is needed, it is strongly recommended that you create a custom PRF
model instead.
If one needs a PRF model from an analytical PSF model, a
more efficient option is to discretize the model on a grid
using :func:`~astropy.convolution.discretize_model` using the
``'oversample'`` or ``'integrate'`` ``mode``. The resulting 2D image
can then be used as the input to ``FittableImageModel`` to create an
ePSF model. This will be *much* faster than using this class.
"""
flux = Parameter(default=1)
x_0 = Parameter(default=0)
y_0 = Parameter(default=0)
def __init__(self, psfmodel, *, renormalize_psf=True, flux=flux.default,
x_0=x_0.default, y_0=y_0.default, xname=None, yname=None,
fluxname=None, **kwargs):
self.psfmodel = psfmodel.copy()
if renormalize_psf:
self._psf_scale_factor = 1.0 / dblquad(self.psfmodel,
-np.inf, np.inf,
lambda x: -np.inf,
lambda x: np.inf)[0]
else:
self._psf_scale_factor = 1
self.xname = xname
self.yname = yname
self.fluxname = fluxname
# these can be used to adjust the integration behavior. Might be
# used in the future to expose how the integration happens
self._dblquadkwargs = {}
super().__init__(n_models=1, x_0=x_0, y_0=y_0, flux=flux, **kwargs)
[docs]
def evaluate(self, x, y, flux, x_0, y_0):
"""
The evaluation function for PRFAdapter.
Parameters
----------
x, y : float or array_like
The coordinates at which to evaluate the model.
flux : float
The total flux of the star.
x_0, y_0 : float
The position of the star.
Returns
-------
evaluated_model : `~numpy.ndarray`
The evaluated model.
"""
if not np.isscalar(flux):
flux = flux[0]
if not np.isscalar(x_0):
x_0 = x_0[0]
if not np.isscalar(y_0):
y_0 = y_0[0]
if self.xname is None:
dx = x - x_0
else:
dx = x
setattr(self.psfmodel, self.xname, x_0)
if self.xname is None:
dy = y - y_0
else:
dy = y
setattr(self.psfmodel, self.yname, y_0)
if self.fluxname is None:
return (flux * self._psf_scale_factor
* self._integrated_psfmodel(dx, dy))
setattr(self.psfmodel, self.yname, flux * self._psf_scale_factor)
return self._integrated_psfmodel(dx, dy)
def _integrated_psfmodel(self, dx, dy):
# infer type/shape from the PSF model. Seems wasteful, but the
# integration step is a *lot* more expensive so its just peanuts
out = np.empty_like(self.psfmodel(dx, dy))
outravel = out.ravel()
for i, (xi, yi) in enumerate(zip(dx.ravel(), dy.ravel(), strict=True)):
outravel[i] = dblquad(self.psfmodel,
xi - 0.5, xi + 0.5,
lambda x: yi - 0.5, lambda x: yi + 0.5,
**self._dblquadkwargs)[0]
return out