Source code for photutils.psf.image_models

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
This module provides image-based PSF models.
"""

import copy
import warnings

import numpy as np
from astropy.modeling import Fittable2DModel, Parameter
from astropy.utils.decorators import deprecated, lazyproperty
from astropy.utils.exceptions import AstropyUserWarning
from scipy.interpolate import RectBivariateSpline

from photutils.aperture import CircularAperture
from photutils.utils._parameters import as_pair

__all__ = ['EPSFModel', 'FittableImageModel', 'ImagePSF']


[docs] class ImagePSF(Fittable2DModel): """ A model for a 2D image PSF. This class takes 2D image data and computes the values of the model at arbitrary locations, including fractional pixel positions, within the image using spline interpolation provided by :py:class:`~scipy.interpolate.RectBivariateSpline`. The model has three model parameters: an image intensity scaling factor (``flux``) which is applied to the input image, and two positional parameters (``x_0`` and ``y_0``) indicating the location of a feature in the coordinate grid on which the model is evaluated. Parameters ---------- data : 2D `~numpy.ndarray` Array containing the 2D image. The length of the x and y axes must both be at least 4. All elements of the input image data must be finite. By default, the PSF peak is assumed to be located at the center of the input image (see the ``origin`` keyword). The array must be normalized so that the total flux of a source is 1.0. This means that the sum of the values in the input image PSF over an infinite grid is 1.0. In practice, the sum of the data values in the input image may be less than 1.0 if the input image only covers a finite region of the PSF. These correction factors can be estimated from the ensquared or encircled energy of the PSF based on the size of the input image. flux : float, optional The total flux of the source, assuming the input image was properly normalized. x_0, y_0 : float The x and y positions of a feature in the image in the output coordinate grid on which the model is evaluated. Typically, this refers to the position of the PSF peak, which is assumed to be located at the center of the input image (see the ``origin`` keyword). origin : tuple of 2 float or None, optional The ``(x, y)`` coordinate with respect to the input image data array that represents the reference pixel of the input data. The reference ``origin`` pixel will be placed at the model ``x_0`` and ``y_0`` coordinates in the output coordinate system on which the model is evaluated. Most typically, the input PSF should be centered in the input image, and thus the origin should be set to the central pixel of the ``data`` array. If the origin is set to `None`, then the origin will be set to the center of the ``data`` array (``(npix - 1) / 2.0``). oversampling : int or array_like (int), optional The integer oversampling factor(s) of the input PSF image. If ``oversampling`` is a scalar then it will be used for both axes. If ``oversampling`` has two elements, they must be in ``(y, x)`` order. fill_value : float, optional The value to use for points outside of the input pixel grid. The default is 0.0. **kwargs : dict, optional Additional optional keyword arguments to be passed to the `astropy.modeling.Model` base class. See Also -------- GriddedPSFModel : A model for a grid of ePSF models. Examples -------- In this simple example, we create a PSF image model from a Circular Gaussian PSF. In this case, one should use the `CircularGaussianPSF` model directly as a PSF model. However, this example demonstrates how to create an image PSF model from an input image. .. plot:: :include-source: import matplotlib.pyplot as plt import numpy as np from photutils.psf import CircularGaussianPSF, ImagePSF gaussian_psf = CircularGaussianPSF(x_0=12, y_0=12, fwhm=3.2) yy, xx = np.mgrid[:25, :25] psf_data = gaussian_psf(xx, yy) psf_model = ImagePSF(psf_data, x_0=12, y_0=12, flux=10) data = psf_model(xx, yy) plt.imshow(data, origin='lower') """ flux = Parameter(default=1, description='Intensity scaling factor of the image.') x_0 = Parameter(default=0, description=('Position of a feature in the image along ' 'the x axis')) y_0 = Parameter(default=0, description=('Position of a feature in the image along ' 'the y axis')) def __init__(self, data, *, flux=flux.default, x_0=x_0.default, y_0=y_0.default, origin=None, oversampling=1, fill_value=0.0, **kwargs): self._validate_data(data) self.data = data self.origin = origin self.oversampling = as_pair('oversampling', oversampling, lower_bound=(0, 1)) self.fill_value = fill_value super().__init__(flux, x_0, y_0, **kwargs) @staticmethod def _validate_data(data): if not isinstance(data, np.ndarray): raise TypeError('Input data must be a 2D numpy array.') if data.ndim != 2: raise ValueError('Input data must be a 2D numpy array.') if not np.all(np.isfinite(data)): raise ValueError('All elements of input data must be finite.') # this is required by RectBivariateSpline for kx=3, ky=3 if np.any(np.array(data.shape) < 4): raise ValueError('The length of the x and y axes must both be at ' 'least 4.') def _cls_info(self): return [('PSF shape (oversampled pixels)', self.data.shape), ('Oversampling', tuple(self.oversampling))] def __str__(self): return self._format_str(keywords=self._cls_info())
[docs] def copy(self): """ Return a copy of this model where only the model parameters are copied. All other copied model attributes are references to the original model. This prevents copying the image data, which may be a large array. This method is useful if one is interested in only changing the model parameters in a model copy. It is used in the PSF photometry classes during model fitting. Use the `deepcopy` method if you want to copy all of the model attributes, including the PSF image data. Returns ------- result : `ImagePSF` A copy of this model with only the model parameters copied. """ newcls = object.__new__(self.__class__) for key, val in self.__dict__.items(): if key in self.param_names: # copy only the parameter values newcls.__dict__[key] = copy.copy(val) else: newcls.__dict__[key] = val return newcls
[docs] def deepcopy(self): """ Return a deep copy of this model. Returns ------- result : `ImagePSF` A deep copy of this model. """ return copy.deepcopy(self)
@property def origin(self): """ A 1D `~numpy.ndarray` (x, y) pixel coordinates within the model's 2D image of the origin of the coordinate system. The reference ``origin`` pixel will be placed at the model ``x_0`` and ``y_0`` coordinates in the output coordinate system on which the model is evaluated. Most typically, the input PSF should be centered in the input image, and thus the origin should be set to the central pixel of the ``data`` array. If the origin is set to `None`, then the origin will be set to the center of the ``data`` array (``(npix - 1) / 2.0``). """ return self._origin @origin.setter def origin(self, origin): if origin is None: origin = (np.array(self.data.shape) - 1.0) / 2.0 origin = origin[::-1] # flip to (x, y) order else: origin = np.asarray(origin) if origin.ndim != 1 or len(origin) != 2: raise ValueError('origin must be 1D and have 2-elements') if not np.all(np.isfinite(origin)): raise ValueError('All elements of origin must be finite') self._origin = origin @lazyproperty def interpolator(self): """ The interpolating spline function. The interpolator is computed with a 3rd-degree `~scipy.interpolate.RectBivariateSpline` (kx=3, ky=3, s=0) using the input image data. The interpolator is used to evaluate the model at arbitrary locations, including fractional pixel positions. Notes ----- This property can be overridden in a subclass to define custom interpolators. """ x = np.arange(self.data.shape[1]) y = np.arange(self.data.shape[0]) # RectBivariateSpline expects the data to be in (x, y) axis order return RectBivariateSpline(x, y, self.data.T, kx=3, ky=3, s=0) def _calc_bounding_box(self): """ Set a bounding box defining the limits of the model. Returns ------- bbox : tuple A bounding box defining the ((y_min, y_max), (x_min, x_max)) limits of the model. """ dy, dx = np.array(self.data.shape) / 2 / self.oversampling # apply the origin shift # if origin is None, the origin is set to the center of the # image and the shift is 0 xshift = np.array(self.data.shape[1] - 1) / 2 - self.origin[0] yshift = np.array(self.data.shape[0] - 1) / 2 - self.origin[1] xshift /= self.oversampling[1] yshift /= self.oversampling[0] return ((self.y_0 - dy + yshift, self.y_0 + dy + yshift), (self.x_0 - dx + xshift, self.x_0 + dx + xshift)) @property def bounding_box(self): """ The bounding box of the model. Examples -------- >>> from photutils.psf import ImagePSF >>> psf_data = np.arange(30, dtype=float).reshape(5, 6) >>> psf_data /= np.sum(psf_data) >>> model = ImagePSF(psf_data, flux=1, x_0=0, y_0=0) >>> model.bounding_box # doctest: +FLOAT_CMP ModelBoundingBox( intervals={ x: Interval(lower=-3.0, upper=3.0) y: Interval(lower=-2.5, upper=2.5) } model=ImagePSF(inputs=('x', 'y')) order='C' ) """ return self._calc_bounding_box()
[docs] def evaluate(self, x, y, flux, x_0, y_0): """ Calculate the value of the image model at the input coordinates for the given model parameters. Parameters ---------- x, y : float or array_like The x and y coordinates at which to evaluate the model. flux : float The total flux of the source, assuming the input image was properly normalized. x_0, y_0 : float The x and y positions of the feature in the image in the output coordinate grid on which the model is evaluated. Returns ------- result : `~numpy.ndarray` The value of the model evaluated at the input coordinates. """ xi = self.oversampling[1] * (np.asarray(x, dtype=float) - x_0) yi = self.oversampling[0] * (np.asarray(y, dtype=float) - y_0) xi += self._origin[0] yi += self._origin[1] evaluated_model = flux * self.interpolator(xi, yi, grid=False) if self.fill_value is not None: # set pixels that are outside the input pixel grid to the # fill_value to avoid extrapolation; these bounds match the # RegularGridInterpolator bounds ny, nx = self.data.shape invalid = (xi < 0) | (xi > nx - 1) | (yi < 0) | (yi > ny - 1) evaluated_model[invalid] = self.fill_value return evaluated_model
[docs] @deprecated('2.0.0', alternative='`ImagePSF`') class FittableImageModel(Fittable2DModel): r""" A fittable image model allowing for intensity scaling and translations. This class takes 2D image data and computes the values of the model at arbitrary locations, including fractional pixel positions, within the image using spline interpolation provided by :py:class:`~scipy.interpolate.RectBivariateSpline`. The fittable model provided by this class has three model parameters: an image intensity scaling factor (``flux``) which is applied to (normalized) image, and two positional parameters (``x_0`` and ``y_0``) indicating the location of a feature in the coordinate grid on which the model is to be evaluated. Parameters ---------- data : 2D `~numpy.ndarray` Array containing the 2D image. flux : float, optional Intensity scaling factor for image data. If ``flux`` is `None`, then the normalization constant will be computed so that the total flux of the model's image data is 1.0. x_0, y_0 : float, optional Position of a feature in the image in the output coordinate grid on which the model is evaluated. normalize : bool, optional Indicates whether or not the model should be build on normalized input image data. If true, then the normalization constant (*N*) is computed so that .. math:: N \cdot C \cdot \sum\limits_{i,j} D_{i,j} = 1, where *N* is the normalization constant, *C* is correction factor given by the parameter ``normalization_correction``, and :math:`D_{i,j}` are the elements of the input image ``data`` array. normalization_correction : float, optional A strictly positive number that represents correction that needs to be applied to model's data normalization (see *C* in the equation in the comments to ``normalize`` for more details). A possible application for this parameter is to account for aperture correction. Assuming model's data represent a PSF to be fitted to some target star, we set ``normalization_correction`` to the aperture correction that needs to be applied to the model. That is, ``normalization_correction`` in this case should be set to the ratio between the total flux of the PSF (including flux outside model's data) to the flux of model's data. Then, best fitted value of the ``flux`` model parameter will represent an aperture-corrected flux of the target star. In the case of aperture correction, ``normalization_correction`` should be a value larger than one, as the total flux, including regions outside of the aperture, should be larger than the flux inside the aperture, and thus the correction is applied as an inversely multiplied factor. origin : tuple, None, optional A reference point in the input image ``data`` array. When origin is `None`, origin will be set at the middle of the image array. If ``origin`` represents the location of a feature (e.g., the position of an intensity peak) in the input ``data``, then model parameters ``x_0`` and ``y_0`` show the location of this peak in an another target image to which this model was fitted. Fundamentally, it is the coordinate in the model's image data that should map to coordinate (``x_0``, ``y_0``) of the output coordinate system on which the model is evaluated. Alternatively, when ``origin`` is set to ``(0, 0)``, then model parameters ``x_0`` and ``y_0`` are shifts by which model's image should be translated in order to match a target image. oversampling : int or array_like (int) The integer oversampling factor(s) of the ePSF relative to the input ``stars`` along each axis. If ``oversampling`` is a scalar then it will be used for both axes. If ``oversampling`` has two elements, they must be in ``(y, x)`` order. fill_value : float, optional The value to be returned by the `evaluate` or ``astropy.modeling.Model.__call__`` methods when evaluation is performed outside the definition domain of the model. **kwargs : dict, optional Additional optional keyword arguments to be passed directly to the `compute_interpolator` method. See `compute_interpolator` for more details. """ flux = Parameter(description='Intensity scaling factor for image data.', default=1.0) x_0 = Parameter(description='X-position of a feature in the image in ' 'the output coordinate grid on which the model is ' 'evaluated.', default=0.0) y_0 = Parameter(description='Y-position of a feature in the image in ' 'the output coordinate grid on which the model is ' 'evaluated.', default=0.0) def __init__(self, data, *, flux=flux.default, x_0=x_0.default, y_0=y_0.default, normalize=False, normalization_correction=1.0, origin=None, oversampling=1, fill_value=0.0, **kwargs): self._fill_value = fill_value self._img_norm = None self._normalization_status = 0 if normalize else 2 self._store_interpolator_kwargs(**kwargs) self._oversampling = as_pair('oversampling', oversampling, lower_bound=(0, 1)) if normalization_correction <= 0: raise ValueError("'normalization_correction' must be strictly " 'positive.') self._normalization_correction = normalization_correction self._normalization_constant = 1.0 / self._normalization_correction self._data = np.array(data, copy=True, dtype=float) if not np.all(np.isfinite(self._data)): raise ValueError("All elements of input 'data' must be finite.") # set input image related parameters: self._ny, self._nx = self._data.shape self._shape = self._data.shape if self._data.size < 1: raise ValueError('Image data array cannot be zero-sized.') # set the origin of the coordinate system in image's pixel grid: self.origin = origin flux = self._initial_norm(flux, normalize) super().__init__(flux, x_0, y_0) # initialize interpolator: self.compute_interpolator(**kwargs) def _initial_norm(self, flux, normalize): if flux is None: if self._img_norm is None: self._img_norm = self._compute_raw_image_norm() flux = self._img_norm self._compute_normalization(normalize) return flux def _compute_raw_image_norm(self): """ Helper function that computes the uncorrected inverse normalization factor of input image data. This quantity is computed as the *sum of all pixel values*. .. note:: This function is intended to be overridden in a subclass if one desires to change the way the normalization factor is computed. """ return np.sum(self._data, dtype=float) def _compute_normalization(self, normalize=True): r""" Helper function that computes (corrected) normalization factor of the original image data. This quantity is computed as the inverse "raw image norm" (or total "flux" of model's image) corrected by the ``normalization_correction``: .. math:: N = 1/(\Phi * C), where :math:`\Phi` is the "total flux" of model's image as computed by `_compute_raw_image_norm` and *C* is the normalization correction factor. :math:`\Phi` is computed only once if it has not been previously computed. Otherwise, the existing (stored) value of :math:`\Phi` is not modified as :py:class:`FittableImageModel` does not allow image data to be modified after the object is created. .. note:: Normally, this function should not be called by the end-user. It is intended to be overridden in a subclass if one desires to change the way the normalization factor is computed. """ self._normalization_constant = 1.0 / self._normalization_correction if normalize: # compute normalization constant so that # N*C*sum(data) = 1: if self._img_norm is None: self._img_norm = self._compute_raw_image_norm() if self._img_norm != 0.0 and np.isfinite(self._img_norm): self._normalization_constant /= self._img_norm self._normalization_status = 0 else: self._normalization_constant = 1.0 self._normalization_status = 1 warnings.warn('Overflow encountered while computing ' 'normalization constant. Normalization ' 'constant will be set to 1.', AstropyUserWarning) else: self._normalization_status = 2 @property def oversampling(self): """ The factor by which the stored image is oversampled. An input to this model is multiplied by this factor to yield the index into the stored image. """ return self._oversampling @property def data(self): """ Get original image data. """ return self._data @property def normalized_data(self): """ Get normalized and/or intensity-corrected image data. """ return self._normalization_constant * self._data @property def normalization_constant(self): """ Get normalization constant. """ return self._normalization_constant @property def normalization_status(self): """ Get normalization status. Possible status values are: * 0: **Performed**. Model has been successfully normalized at user's request. * 1: **Failed**. Attempt to normalize has failed. * 2: **NotRequested**. User did not request model to be normalized. """ return self._normalization_status @property def normalization_correction(self): """ Set/Get flux correction factor. .. note:: When setting correction factor, model's flux will be adjusted accordingly such that if this model was a good fit to some target image before, then it will remain a good fit after correction factor change. """ return self._normalization_correction @normalization_correction.setter def normalization_correction(self, normalization_correction): old_cf = self._normalization_correction self._normalization_correction = normalization_correction self._compute_normalization(normalize=self._normalization_status != 2) # adjust model's flux so that if this model was a good fit to # some target image, then it will remain a good fit after # correction factor change: self.flux *= normalization_correction / old_cf @property def shape(self): """ A tuple of dimensions of the data array in numpy style (ny, nx). """ return self._shape @property def nx(self): """ Number of columns in the data array. """ return self._nx @property def ny(self): """ Number of rows in the data array. """ return self._ny @property def origin(self): """ A tuple of ``x`` and ``y`` coordinates of the origin of the coordinate system in terms of pixels of model's image. When setting the coordinate system origin, a tuple of two integers or floats may be used. If origin is set to `None`, the origin of the coordinate system will be set to the middle of the data array (``(npix-1)/2.0``). .. warning:: Modifying ``origin`` will not adjust (modify) model's parameters ``x_0`` and ``y_0``. """ return (self._x_origin, self._y_origin) @origin.setter def origin(self, origin): if origin is None: self._x_origin = (self._nx - 1) / 2.0 self._y_origin = (self._ny - 1) / 2.0 elif hasattr(origin, '__iter__') and len(origin) == 2: self._x_origin, self._y_origin = origin else: raise TypeError('Parameter "origin" must be either None or an ' 'iterable with two elements.') @property def x_origin(self): """ X-coordinate of the origin of the coordinate system. """ return self._x_origin @property def y_origin(self): """ Y-coordinate of the origin of the coordinate system. """ return self._y_origin @property def fill_value(self): """ Fill value to be returned for coordinates outside of the domain of definition of the interpolator. If ``fill_value`` is `None`, then values outside of the domain of definition are the ones returned by the interpolator. """ return self._fill_value @fill_value.setter def fill_value(self, fill_value): self._fill_value = fill_value def _store_interpolator_kwargs(self, **kwargs): """ Store interpolator keyword arguments. This function should be called in a subclass whenever model's interpolator is (re)computed. """ self._interpolator_kwargs = copy.deepcopy(kwargs) @property def interpolator_kwargs(self): """ Get current interpolator's arguments used when interpolator was created. """ return self._interpolator_kwargs
[docs] def compute_interpolator(self, **kwargs): """ Compute/define the interpolating spline. This function can be overridden in a subclass to define custom interpolators. Parameters ---------- **kwargs : dict, optional Additional optional keyword arguments: * **degree** : int, tuple, optional Degree of the interpolating spline. A tuple can be used to provide different degrees for the X- and Y-axes. Default value is degree=3. * **s** : float, optional Non-negative smoothing factor. Default value s=0 corresponds to interpolation. See :py:class:`~scipy.interpolate.RectBivariateSpline` for more details. Notes ----- * When subclassing :py:class:`FittableImageModel` for the purpose of overriding :py:func:`compute_interpolator`, the :py:func:`evaluate` may need to overridden as well depending on the behavior of the new interpolator. In addition, for improved future compatibility, make sure that the overriding method stores keyword arguments ``kwargs`` by calling ``_store_interpolator_kwargs`` method. * Use caution when modifying interpolator's degree or smoothness in a computationally intensive part of the code as it may decrease code performance due to the need to recompute interpolator. """ if 'degree' in kwargs: degree = kwargs['degree'] if hasattr(degree, '__iter__') and len(degree) == 2: degx = int(degree[0]) degy = int(degree[1]) else: degx = int(degree) degy = int(degree) if degx < 0 or degy < 0: raise ValueError('Interpolator degree must be a non-negative ' 'integer') else: degx = 3 degy = 3 smoothness = kwargs.get('s', 0) x = np.arange(self._nx, dtype=float) y = np.arange(self._ny, dtype=float) self.interpolator = RectBivariateSpline( x, y, self._data.T, kx=degx, ky=degy, s=smoothness ) self._store_interpolator_kwargs(**kwargs)
[docs] def evaluate(self, x, y, flux, x_0, y_0, *, use_oversampling=True): """ Calculate the value of the image model at the input coordinates for the given model parameters. Parameters ---------- x, y : float or array_like The x and y coordinates at which to evaluate the model. flux : float The total flux of the source. x_0, y_0 : float The x and y positions of the feature in the image in the output coordinate grid on which the model is evaluated. use_oversampling : bool, optional Whether to use the oversampling factor to calculate the model pixel indices. The default is `True`, which means the input indices will be multiplied by this factor. Returns ------- evaluated_model : `~numpy.ndarray` The evaluated model. """ if use_oversampling: xi = self._oversampling[1] * (np.asarray(x) - x_0) yi = self._oversampling[0] * (np.asarray(y) - y_0) else: xi = np.asarray(x) - x_0 yi = np.asarray(y) - y_0 xi = xi.astype(float) yi = yi.astype(float) xi += self._x_origin yi += self._y_origin f = flux * self._normalization_constant evaluated_model = f * self.interpolator.ev(xi, yi) if self._fill_value is not None: # find indices of pixels that are outside the input pixel grid and # set these pixels to the 'fill_value': invalid = (((xi < 0) | (xi > self._nx - 1)) | ((yi < 0) | (yi > self._ny - 1))) evaluated_model[invalid] = self._fill_value return evaluated_model
class _LegacyEPSFModel(Fittable2DModel): """ This class will be removed when the deprecated EPSFModel is removed, which will require the EPSFBuilder class to be rewritten/refactored/replaced. A class that models an effective PSF (ePSF). The EPSFModel is normalized such that the sum of the PSF over the (undersampled) pixels within the input ``norm_radius`` is 1.0. This means that when the EPSF is fit to stars, the resulting flux corresponds to aperture photometry within a circular aperture of radius ``norm_radius``. While this class is a subclass of `FittableImageModel`, it is very similar. The primary differences/motivation are a few additional parameters necessary specifically for ePSFs. Parameters ---------- data : 2D `~numpy.ndarray` Array containing the 2D image. flux : float, optional Intensity scaling factor for image data. x_0, y_0 : float, optional Position of a feature in the image in the output coordinate grid on which the model is evaluated. normalize : bool, optional Indicates whether or not the model should be build on normalized input image data. normalization_correction : float, optional A strictly positive number that represents correction that needs to be applied to model's data normalization. origin : tuple, None, optional A reference point in the input image ``data`` array. When origin is `None`, origin will be set at the middle of the image array. oversampling : int or array_like (int) The integer oversampling factor(s) of the ePSF relative to the input ``stars`` along each axis. If ``oversampling`` is a scalar then it will be used for both axes. If ``oversampling`` has two elements, they must be in ``(y, x)`` order. fill_value : float, optional The value to be returned when evaluation is performed outside the domain of the model. norm_radius : float, optional The radius inside which the ePSF is normalized by the sum over undersampled integer pixel values inside a circular aperture. **kwargs : dict, optional Additional optional keyword arguments to be passed directly to the `compute_interpolator` method. See `compute_interpolator` for more details. """ flux = Parameter(description='Intensity scaling factor for image data.', default=1.0) x_0 = Parameter(description='X-position of a feature in the image in ' 'the output coordinate grid on which the model is ' 'evaluated.', default=0.0) y_0 = Parameter(description='Y-position of a feature in the image in ' 'the output coordinate grid on which the model is ' 'evaluated.', default=0.0) def __init__(self, data, *, flux=flux.default, x_0=x_0.default, y_0=y_0.default, normalize=False, normalization_correction=1.0, origin=None, oversampling=1, fill_value=0.0, norm_radius=5.5, **kwargs): self._norm_radius = norm_radius self._fill_value = fill_value self._img_norm = None self._normalization_status = 0 if normalize else 2 self._store_interpolator_kwargs(**kwargs) self._oversampling = as_pair('oversampling', oversampling, lower_bound=(0, 1)) if normalization_correction <= 0: raise ValueError("'normalization_correction' must be strictly " 'positive.') self._normalization_correction = normalization_correction self._normalization_constant = 1.0 / self._normalization_correction self._data = np.array(data, copy=True, dtype=float) if not np.all(np.isfinite(self._data)): raise ValueError("All elements of input 'data' must be finite.") # set input image related parameters: self._ny, self._nx = self._data.shape self._shape = self._data.shape if self._data.size < 1: raise ValueError('Image data array cannot be zero-sized.') # set the origin of the coordinate system in image's pixel grid: self.origin = origin flux = self._initial_norm(flux, normalize) super().__init__(flux, x_0, y_0) # initialize interpolator: self.compute_interpolator(**kwargs) def _initial_norm(self, flux, normalize): if flux is None: if self._img_norm is None: self._img_norm = self._compute_raw_image_norm() flux = self._img_norm if normalize: self._compute_normalization() else: self._img_norm = self._compute_raw_image_norm() return flux def _compute_raw_image_norm(self): """ Compute the normalization of input image data as the flux within a given radius. """ xypos = (self._nx / 2.0, self._ny / 2.0) # TODO: generalize "radius" (ellipse?) is oversampling is # different along x/y axes radius = self._norm_radius * self.oversampling[0] aper = CircularAperture(xypos, r=radius) flux, _ = aper.do_photometry(self._data, method='exact') return flux[0] / np.prod(self.oversampling) def _compute_normalization(self, normalize=True): """ Helper function that computes (corrected) normalization factor of the original image data. For the ePSF this is defined as the sum over the inner N (default=5.5) pixels of the non-oversampled image. Will renormalize the data to the value calculated. """ if normalize: if self._img_norm is None: if np.sum(self._data) == 0: self._img_norm = 1 else: self._img_norm = self._compute_raw_image_norm() if self._img_norm != 0.0 and np.isfinite(self._img_norm): self._data /= (self._img_norm * self._normalization_correction) self._normalization_status = 0 else: self._normalization_status = 1 self._img_norm = 1 warnings.warn('Overflow encountered while computing ' 'normalization constant. Normalization ' 'constant will be set to 1.', AstropyUserWarning) else: self._normalization_status = 2 @property def normalized_data(self): """ Overloaded dummy function that also returns self._data, as the normalization occurs within _compute_normalization in EPSFModel, and as such self._data will sum, accounting for under/oversampled pixels, to 1/self._normalization_correction. """ return self._data @property def oversampling(self): """ The factor by which the stored image is oversampled. An input to this model is multiplied by this factor to yield the index into the stored image. """ return self._oversampling @property def data(self): """ Get original image data. """ return self._data @property def normalization_constant(self): """ Get normalization constant. """ return self._normalization_constant @property def normalization_status(self): """ Get normalization status. Possible status values are: * 0: **Performed**. Model has been successfully normalized at user's request. * 1: **Failed**. Attempt to normalize has failed. * 2: **NotRequested**. User did not request model to be normalized. """ return self._normalization_status @property def normalization_correction(self): """ Set/Get flux correction factor. .. note:: When setting correction factor, model's flux will be adjusted accordingly such that if this model was a good fit to some target image before, then it will remain a good fit after correction factor change. """ return self._normalization_correction @normalization_correction.setter def normalization_correction(self, normalization_correction): old_cf = self._normalization_correction self._normalization_correction = normalization_correction self._compute_normalization(normalize=self._normalization_status != 2) # adjust model's flux so that if this model was a good fit to # some target image, then it will remain a good fit after # correction factor change: self.flux *= normalization_correction / old_cf @property def shape(self): """ A tuple of dimensions of the data array in numpy style (ny, nx). """ return self._shape @property def nx(self): """ Number of columns in the data array. """ return self._nx @property def ny(self): """ Number of rows in the data array. """ return self._ny @property def origin(self): """ A tuple of ``x`` and ``y`` coordinates of the origin of the coordinate system in terms of pixels of model's image. When setting the coordinate system origin, a tuple of two integers or floats may be used. If origin is set to `None`, the origin of the coordinate system will be set to the middle of the data array (``(npix-1)/2.0``). .. warning:: Modifying ``origin`` will not adjust (modify) model's parameters ``x_0`` and ``y_0``. """ return (self._x_origin, self._y_origin) @origin.setter def origin(self, origin): if origin is None: self._x_origin = (self._nx - 1) / 2.0 / self.oversampling[1] self._y_origin = (self._ny - 1) / 2.0 / self.oversampling[0] elif (hasattr(origin, '__iter__') and len(origin) == 2): self._x_origin, self._y_origin = origin else: raise TypeError('Parameter "origin" must be either None or an ' 'iterable with two elements.') @property def x_origin(self): """ X-coordinate of the origin of the coordinate system. """ return self._x_origin @property def y_origin(self): """ Y-coordinate of the origin of the coordinate system. """ return self._y_origin @property def fill_value(self): """ Fill value to be returned for coordinates outside of the domain of definition of the interpolator. If ``fill_value`` is `None`, then values outside of the domain of definition are the ones returned by the interpolator. """ return self._fill_value @fill_value.setter def fill_value(self, fill_value): self._fill_value = fill_value def _store_interpolator_kwargs(self, **kwargs): """ Store interpolator keyword arguments. This function should be called in a subclass whenever model's interpolator is (re)computed. """ self._interpolator_kwargs = copy.deepcopy(kwargs) @property def interpolator_kwargs(self): """ Get current interpolator's arguments used when interpolator was created. """ return self._interpolator_kwargs def compute_interpolator(self, **kwargs): """ Compute/define the interpolating spline. This function can be overridden in a subclass to define custom interpolators. Parameters ---------- **kwargs : dict, optional Additional optional keyword arguments: * **degree** : int, tuple, optional Degree of the interpolating spline. A tuple can be used to provide different degrees for the X- and Y-axes. Default value is degree=3. * **s** : float, optional Non-negative smoothing factor. Default value s=0 corresponds to interpolation. See :py:class:`~scipy.interpolate.RectBivariateSpline` for more details. Notes ----- * When subclassing :py:class:`FittableImageModel` for the purpose of overriding :py:func:`compute_interpolator`, the :py:func:`evaluate` may need to overridden as well depending on the behavior of the new interpolator. In addition, for improved future compatibility, make sure that the overriding method stores keyword arguments ``kwargs`` by calling ``_store_interpolator_kwargs`` method. * Use caution when modifying interpolator's degree or smoothness in a computationally intensive part of the code as it may decrease code performance due to the need to recompute interpolator. """ if 'degree' in kwargs: degree = kwargs['degree'] if hasattr(degree, '__iter__') and len(degree) == 2: degx = int(degree[0]) degy = int(degree[1]) else: degx = int(degree) degy = int(degree) if degx < 0 or degy < 0: raise ValueError('Interpolator degree must be a non-negative ' 'integer') else: degx = 3 degy = 3 smoothness = kwargs.get('s', 0) # Interpolator must be set to interpolate on the undersampled # pixel grid, going from 0 to len(undersampled_grid) x = np.arange(self._nx, dtype=float) / self.oversampling[1] y = np.arange(self._ny, dtype=float) / self.oversampling[0] self.interpolator = RectBivariateSpline( x, y, self._data.T, kx=degx, ky=degy, s=smoothness ) self._store_interpolator_kwargs(**kwargs) def evaluate(self, x, y, flux, x_0, y_0): """ Calculate the value of the image model at the input coordinates for the given model parameters. Parameters ---------- x, y : float or array_like The x and y coordinates at which to evaluate the model. flux : float The total flux of the source. x_0, y_0 : float The x and y positions of the feature in the image in the output coordinate grid on which the model is evaluated. Returns ------- evaluated_model : `~numpy.ndarray` The evaluated model. """ xi = np.asarray(x) - x_0 + self._x_origin yi = np.asarray(y) - y_0 + self._y_origin evaluated_model = flux * self.interpolator.ev(xi, yi) if self._fill_value is not None: # find indices of pixels that are outside the input pixel # grid and set these pixels to the 'fill_value': invalid = (((xi < 0) | (xi > (self._nx - 1) / self.oversampling[1])) | ((yi < 0) | (yi > (self._ny - 1) / self.oversampling[0]))) evaluated_model[invalid] = self._fill_value return evaluated_model
[docs] @deprecated('2.0.0', alternative='`ImagePSF`') class EPSFModel(_LegacyEPSFModel): """ A class that models an effective PSF (ePSF). The EPSFModel is normalized such that the sum of the PSF over the (undersampled) pixels within the input ``norm_radius`` is 1.0. This means that when the EPSF is fit to stars, the resulting flux corresponds to aperture photometry within a circular aperture of radius ``norm_radius``. While this class is a subclass of `FittableImageModel`, it is very similar. The primary differences/motivation are a few additional parameters necessary specifically for ePSFs. Parameters ---------- data : 2D `~numpy.ndarray` Array containing the 2D image. flux : float, optional Intensity scaling factor for image data. x_0, y_0 : float, optional Position of a feature in the image in the output coordinate grid on which the model is evaluated. normalize : bool, optional Indicates whether or not the model should be build on normalized input image data. normalization_correction : float, optional A strictly positive number that represents correction that needs to be applied to model's data normalization. origin : tuple, None, optional A reference point in the input image ``data`` array. When origin is `None`, origin will be set at the middle of the image array. oversampling : int or array_like (int) The integer oversampling factor(s) of the ePSF relative to the input ``stars`` along each axis. If ``oversampling`` is a scalar then it will be used for both axes. If ``oversampling`` has two elements, they must be in ``(y, x)`` order. fill_value : float, optional The value to be returned when evaluation is performed outside the domain of the model. norm_radius : float, optional The radius inside which the ePSF is normalized by the sum over undersampled integer pixel values inside a circular aperture. **kwargs : dict, optional Additional optional keyword arguments to be passed directly to the "compute_interpolator" method. See "compute_interpolator" for more details. """ flux = Parameter(description='Intensity scaling factor for image data.', default=1.0) x_0 = Parameter(description='X-position of a feature in the image in ' 'the output coordinate grid on which the model is ' 'evaluated.', default=0.0) y_0 = Parameter(description='Y-position of a feature in the image in ' 'the output coordinate grid on which the model is ' 'evaluated.', default=0.0) def __init__(self, data, *, flux=flux.default, x_0=x_0.default, y_0=y_0.default, normalize=True, normalization_correction=1.0, origin=None, oversampling=1, fill_value=0.0, norm_radius=5.5, **kwargs): super().__init__(data=data, flux=flux, x_0=x_0, y_0=y_0, normalize=normalize, normalization_correction=normalization_correction, origin=origin, oversampling=oversampling, fill_value=fill_value, norm_radius=norm_radius, **kwargs)