pyimcom.psfutil
Utilities for PSFs and system matrices.
Classes
- OutPSF
Simple target output PSF models.
- PSFGrp
Group of PSFs attached to either an InStamp or a Block.
- PSFOvl
Overlap between two PSFGrp instances or a PSFGrp instance and itself.
- SysMatA
System matrix A attached to a coadd.Block instance.
- SysMatB
System matrix B attached to a coadd.Block instance.
Attributes
Classes
Selector for PSF interpolations. |
|
Simple target output PSF models (for testing or outputs). |
|
Group of PSFs. |
|
Overlap between two PSFGrp instances or a PSFGrp instance and itself. |
|
System matrix A attached to a coadd.Block instance. |
|
System matrix B attached to a coadd.Block instance. |
Module Contents
- class OutPSF[source]
Simple target output PSF models (for testing or outputs).
- psf_cplx_airy()[source]
Airy function with extra features printed on top to make it messier (staticmethod).
- psf_get_fwhm()
Compute full width at half maximum of a PSF (staticmethod).
- psf_get_inv_width()
Compute shear-invariant width of a PSF (staticmethod).
- static psf_gaussian(n: int, sigmax: float, sigmay: float) numpy.array[source]
Gaussian spot, n x n, given sigma, centered (useful for testing).
- Parameters:
n (int) – Size of output to generate.
sigmax (float) – 1 sigma width in the x-direction.
sigmay (float) – 1 sigma width in the y-direction.
- Returns:
Output image, shape (n, n).
- Return type:
np.array of flat
- static psf_simple_airy(n: int, ldp: float, obsc: float = 0.0, tophat_conv: float = 0.0, sigma: float = 0.0) numpy.array[source]
Airy spot, optionally obscured and convolved with a Gaussian and tophat.
Airy spot, n x n, with lambda/D = ldp pixels, and convolved with a tophat (square, full width tophat_conv) and Gaussian (sigma) and linear obscuration factor (obsc).
- Parameters:
n (int) – Size of output to generate.
ldp (float) – Diffraction spot width lambda/D in pixels.
obsc (float, optional) – Fractional linear obscuration.
tophat_conv (float, optional) – Convolution with a square top-hat of the given full width.
sigma (float, optional) – Convolution with a Gaussian of the given 1 sigma width.
- Returns:
The output spot. Shape (n, n).
- Return type:
np.array
Notes
The result is centered on (n-1)/2, (n-1)/2 (so on a pixel if n is odd and a corner if n is even).
The image is normalized to sum to unity if analytically extended to infinity (so with a stamp enclosing 90% of the energy, will sum to 0.9).
- static psf_cplx_airy(n: int, ldp: float, tophat_conv: float = 0.0, sigma: float = 0.0, features: int = 0) numpy.array[source]
somewhat messier Airy function with a few diffraction features printed on ‘features’ is an integer that can be added. everything is band limited
- Parameters:
n (int) – Size of output to generate.
ldp (float) – Diffraction spot width lambda/D in pixels.
tophat_conv (float, optional) – Convolution with a square top-hat of the given full width.
sigma (float, optional) – Convolution with a Gaussian of the given 1 sigma width.
features (int, optional) – Flags controlling which messy features to add.
- Returns:
The output spot. Shape (n, n).
- Return type:
np.array
See also
psf_simple_airyVersion without messy added features.
- static iD5512C_getw(w: numpy.array, fh: float) None[source]
No-Numba version of routine.iD5512C_getw.
The Numba version may not work in a Jupyter environment.
- Parameters:
w (np.array, shape : (10,)) – Interpolation weights in one direction.
fh (float) – ‘xfh’ and ‘yfh’ with 1/2 subtracted.
- Return type:
None.
- static get_psf_fwhm(psf: numpy.array, visualize: bool = False) float[source]
Compute FWHM of a PSF.
This function assumes that the given PSF is azimuthally symmetric.
- Parameters:
psf (np.array, shape : (ny, nx)) – PSF of which the FWHM is to be computed. Typically something returned by PSFGrp._get_outpsf.
visualize (bool, optional) – Whether to visualize PSF and FWHM. The default is False.
- Returns:
FWHM of given PSF in units of pixels.
- Return type:
float
- static get_psf_inv_width(psf: numpy.array) float[source]
Compute shear-invariant width of a PSF.
- Parameters:
psf (np.array, shape : (ny, nx)) – PSF of which the shear-invariant width is to be computed. Typically something returned by PSFGrp._get_outpsf.
- Returns:
Shear-invariant width of given PSF in units of pixels.
- Return type:
float
- class PSFGrp(in_or_out: bool = True, inst=None, blk=None, verbose: bool = False, visualize: bool = False)[source]
Group of PSFs.
Either a group of input PSFs attached to an coadd.InStamp instance or a group of output PSFs attached to a coadd.Block instance.
- Parameters:
in_or_out (bool, optional) – True if input PSF group, False if output PSF group.
inst (pyimcom.coadd.InStamp) – Input stamp, must be provided when in_or_out is True.
blk (pyimcom.coadd.Block) – Block, must be provided when in_or_out is False.
verbose (bool, optional) – Whether to print additional information.
visualize (bool, optional) – Whether to visualize PSFs and sampling positions.
- classmethod setup(npixpsf: int = 48, oversamp: int = 8, dtheta: float = 0.025 / 3600, psfsplit: bool = False) None[source]
Set up class attributes.
- Parameters:
npixpsf (int, optional) – Size of PSF postage stamp in native pixels.
oversamp (int, optional) – PSF oversampling factor relative to native pixel scale.
dtheta (float, optional) – Output pixel scale in degrees. The default is 0.025/3600 (corresponds to 0.025 arcsec)).
psfsplit (bool, optional) – Is PSF splitting implemented in this run?
- Return type:
None
- static visualize_psf(psf_: numpy.array, yxco: numpy.array, xctr: float, yctr: float) None[source]
Visualize a single PSF, together with sampling positions.
- Parameters:
psf (np.array) – A single PSF. Must be 2-dimensional.
yxco (np.array) – Sampling positions, with (0, 0) at the PSF center. Must be 3-dimensional. `yxco`[0,:,:] represents x-positions and `yxco`[1,:,:] represents y-positions.
xctr (float) – Position of the PSF center.
yctr (float) – Position of the PSF center.
- Return type:
None
- _sample_psf(idx, psf, outpix2world2inpix=None, visualize=False)[source]
Perform interpolations to sample a single PSF or a set of PSFs.
- Parameters:
idx (int or None) – Index of the single PSF (int). If None, sample a set of PSFs at the same time.
psf (np.array) – PSF(s) to be sampled. Shape (n_psf, ny, nx).
outpix2world2inpix (method, optional) – InImage outpix2world2inpix of the appropriate InImage instance. The default is None, meaning not to rotate the sampling positions.
visualize (bool, optional) – Whether to visualize PSFs and sampling positions.
- Return type:
None
- _build_inpsfgrp(visualize: bool = False) None[source]
Build a group of input PSFs.
- Parameters:
visualize (bool, optional) – Whether to visualize PSFs and sampling positions.
- Return type:
None
- static _get_outpsf(outpsf: str = 'AIRYOBSC', extrasmooth: float = 0.0, use_filter: int = 4) numpy.array[source]
Get an output PSF specified by configuration.
- Parameters:
outpsf (str, optional) – Target output PSF type. Options are ‘GAUSSIAN’, ‘AIRYOBSC’, and ‘AIRYUNOBSC’.
extrasmooth (float, optional) – Target output PSF extra smearing. The default is 0.0.
use_filter (int, optional) – Which filter to use (0..10).
- Returns:
np.array, shape – Output PSF specified by configuration.
- Return type:
(PSFGrp.nsamp+1, PSFGrp.nsamp+1)
- _build_outpsfgrp(visualize: bool = False) None[source]
Build a group of output PSFs.
What output PSFs to include is set by the configuration file.
- Parameters:
visualize (bool, optional) – Whether to visualize PSFs and sampling positions.
- Return type:
None
- static accel_pad_and_rfft2(psf_arr: numpy.array) numpy.array[source]
Accelerated version of zero-padding and rfft2.
For nsamp = 537 and nfft = 1280, this saves about 20% of the time compared to simply using rfft2.
- Parameters:
psf_arr (np.array) – PSF array to be padded and Fourier transformed. The shape should be (…, PSFGrp.nsamp, PSFGrp.nsamp).
- Returns:
Real Fourier transform results of psf_arr. The shape is (…, PSFGrp.nfft, PSFGrp.nfft//2+1).
- Return type:
np.array
Notes
Diagram of the steps involved here:
# +---+ +---+ # | | fft |***| # +---+ pad +-------+ rfft +---+ pad | | ===> |***| # |***| ===> |*** | ===> |***| ===> |***| |***| # +---+ +-------+ +---+ +---+ +---+ # psf_arr pad_m1 pad_m2 res
Note: m1 stands for axis=-1, and m2 stands for axis=-2.
- class PSFOvl(psfgrp1: PSFGrp, psfgrp2: PSFGrp = None, verbose: bool = False, visualize: bool = False)[source]
Overlap between two PSFGrp instances or a PSFGrp instance and itself.
- Parameters:
psfgrp1 (pyimcom.psfutil.PSFGrp) – The first PSFGrp instance.
psfgrp2 (pyimcom.psfutil.PSFGrp or None, optional) – The second PSFGrp instance. The default is None, indicating the self-overlap of psfgrp1.
verbose (bool, optional) – Whether to print additional information.
visualize (bool, optional) – Whether to visualize the PSF overlap array.
- classmethod setup(flat_penalty: float = 1e-07) None[source]
Set up class attribute.
This must be run after
PSFGrp.setup.- Parameters:
flat_penalty (float, optional) – Amount by which to penalize having different contributions to the output from different input images.
- Return type:
None
- _idx_square2triangle(idx1: int, idx2: int) int[source]
Convert a 2D square index to a 1D triangle index.
- Parameters:
idx1 (int, int) – 2D square index.
idx2 (int, int) – 2D square index.
- Returns:
1D triangle index.
- Return type:
int
Notes
The motivation is that, in the input self-overlap case, the overlap between psf_arr[idx2] and psf_arr[idx1] is simply the flipped version of the overlap between psf_arr[idx1] and psf_arr[idx2]. Therefore, we only need to compute the latter, and the results are stored in a np.array of shape (n_psf*(n_psf+1)//2, nsamp, nsamp).
For example, in the n_psf == 3 case:
# \idx2 0 1 2 # idx1 +---+---+---+ # 0 | 0 | 1 | 2 | # +---+---+---+ # 1 | | 3 | 4 | # +---+---+---+ # 2 | | | 5 | # +---+---+---+
- static accel_irfft2_and_extract(ovl_rft: numpy.array) numpy.array[source]
Accelerated version of irfft2 and extraction (used in FFT-based convolution).
- Parameters:
ovl_rft (np.array) – Real Fourier transform of the PSF overlap array we want. The shape is (…, PSFGrp.nfft, PSFGrp.nfft//2+1).
- Returns:
The PSF overlap array we want. The shape is (…, nsamp, nsamp)
- Return type:
np.array
Notes
This acceleration is based on the fact that we only need less than (1/2) x (1/2) of the irfft2 results. For nsamp = 537 and nfft = 1280, this saves about 40% of the time compared to simply using irfft2 and ifftshift.
Because of the ifftshift part, this function is not the inverse function of PSFGrp.accel_pad_and_rfft2.
Diagram of the steps involved here (ext. stands for extraction):
# +---+ +---+ # |***| ifft |***| # |***| ===> |***| ext. +---+ irfft+-------+ ext. +---+ # |***| |***| ===> |***| ===> |*******| ===> |***| # +---+ +---+ +---+ +-------+ +---+ # ovl_rft ift_m2 ovl_m2 ift_m1 ovl_m1
Note: m2 stands for axis=-2, and m1 stands for axis=-1.
If PSFOvl.nsamp is at least half of PSFGrp.nfft, then reverts to using irfft2 and ifftshift.
- _build_psfovl(visualize: bool = False) None[source]
Build the PSF overlap array.
- Parameters:
visualize (bool, optional) – Whether to visualize the PSF overlap array.
- Return type:
None
- __call__(st1, st2=None, visualize=False) numpy.array[source]
Wrapper for the C interpolators.
- Parameters:
st1 (pyimcom.coadd.InStamp) – First input stamp.
st2 (pyimcom.coadd.InStamp or coadd.OutStamp or None) – Second input stamp. The default is None, indicating diagonal blocks of the A matrix (2nd stamp same as 1st).
visualize (bool, optional) – Whether to visualize overlap arrays and sampling positions.
- Returns:
The shape depends on the nature of this PSFOvl instance and the input.
- Return type:
np.array
- _call_ii_cross(st1, st2, visualize=False) numpy.array[source]
Interpolations in the input-input cross-overlap case.
- Parameters:
st1 (pyimcom.coadd.InStamp) – First input stamp.
st2 (pyimcom.coadd.InStamp) – Second input stamp.
visualize (bool, optional) – Whether to visualize overlap arrays and sampling positions.
- Returns:
A 2D array of the cross-correlation between the two input stamps. The shape is (st1.pix_cumsum[-1], st2.pix_cumsum[-1]).
- Return type:
np.array
- _call_io_cross(st1, st2, visualize=False)[source]
Interpolations in the input-output cross-overlap case.
- Parameters:
st1 (pyimcom.coadd.InStamp) – Input stamp.
st2 (pyimcom.coadd.OutStamp) – Output stamp.
visualize (bool, optional) – Whether to visualize overlap arrays and sampling positions.
- Returns:
A 3D image of the cross-correlation of the input and output PSFs, with axis=0 indicating which output PSF is used. The shape is either (self.grp2.n_psf, n_outpix, st1.pix_cumsum[-1]) or (self.grp2.n_psf, n_outpix, selection.shape[0]).
- Return type:
np.array
- _call_ii_self(st1, st2, visualize=False)[source]
Interpolations in the input-input self-overlap case.
Note that the self-overlap of an input PSF group can be used to compute either the block-diagonal submatrices for a single input stamp or the submatrices for a pair of input stamps within a 2x2 group.
- Parameters:
st1 (pyimcom.coadd.InStamp) – First input stamp.
st2 (pyimcom.coadd.InStamp or None, optional) – Second input stamp. The default is None, indicating diagonal blocks of the A matrix.
visualize (bool, optional) – Whether to visualize overlap arrays and sampling positions.
- Returns:
A 2D image of the cross-correlation. The shape is (st1.pix_cumsum[-1], st2.pix_cumsum[-1]).
- Return type:
np.array
- class SysMatA(blk)[source]
System matrix A attached to a coadd.Block instance.
The symtem matrix A is defined in Rowe+ 2011 Equation (17).
- Parameters:
blk (pyimcom.coadd.Block) – The Block instance to which this SysMatA instance is attached.
- static ji_st2psf(ji_st: int, int)[source]
Convert InStamp index to PSFGrp index.
More precisely, convert index of the InStamp of interest to that of its neighbor to which the PSFGrp in that 2x2 group is attached. For example, any element of {(0, 0), (0, 1), (1, 0), (1, 1)} -> (0, 0).
- Parameters:
ji_st ((int, int)) – Index of the InStamp of interest.
- Returns:
Index of the InStamp which harbors the PSFGrp.
- Return type:
(int, int)
- static shift_ji_st(ji_st: int, int, dji_st: int, int)[source]
Tool function for tuple addition.
For our purposes, this function shifts an InStamp index in each direction by either 0 or 1.
- Parameters:
ji_st ((int, int)) – InStamp index to be shifted.
dji_st ((int, int)) – Shift in InStamp index.
- Returns:
Shifted InStamp index.
- Return type:
(int, int)
- static iisubmat_dist(ji_st1: int, int, ji_st2: int, int)[source]
Calculate the index for iisubmats_ref.
Specifically, calculate the “distance” between two input postage stamps in order to store the reference counts in the 3D array iisubmats_ref using indices (`ji_st1`[0], `ji_st1`[1], dist), which this static method returns.
- Parameters:
ji_st1 ((int, int)) – Index of the first InStamp.
ji_st2 ((int, int)) – Index of the second InStamp.
- Returns:
Index of the first InStamp combined with the “distance,” which can be directly used as the index of iisubmats_ref. Returns None when the “distance” is out of range.
- Return type:
(int, int, int) or None
Notes
The “distance” is defined as follows:
# +---+---+---+---+---+ # | 8 | 9 | 10| 11| 12| # +---+---+---+---+---+ # | 3 | 4 | 5 | 6 | 7 | # +---+---+---+---+---+ # | | | 0*| 1 | 2 | # +---+---+---+---+---+
where “*” denotes ji_st1. Note that ji_st1 must precede ji_st2. If the “distance” is out of range, this function returns None.
- _compute_iisubmats(ji_st1: int, int, ji_st2: int, int, sim_mode: bool = False, verbose: bool = False) None[source]
Make input-input PSFOvl and compute A submatrices.
This method is only called in SysMatA.get_iisubmat. When an A submatrix is requested from an OutStamp but not in iisubmats, get_iisubmat calls this method, which makes the input-input PSFOvl, computes all the needed A submatrices which are supposed to be computed using that particular PSFOvl instance, and clears the PSFOvl.
Note that removing the huge PSF overlap arrays is important! Otherwise pyimcom would demand much more memory.
To determine whether a specific system submatrix is needed, a Block instance first runs the postage coaddition framework in sim_mode to count references to all possible system submatrices.
- Parameters:
ji_st1 ((int, int)) – Index of the first InStamp.
ji_st2 ((int, int)) – Index of the second InStamp.
sim_mode (bool, optional) – Whether to count references without actually computing submatrices.
verbose (bool, optional) – Whether to print additional information.
- Return type:
None
- get_iisubmat(ji_st1: int, int, ji_st2: int, int, sim_mode: bool = False, ji_st_out: int, int = None) numpy.array[source]
Return the requested A submatrix.
This is the only public interface of this class. For the sim_mode, see the docstring of SysMatA._compute_iisubmats.
- Parameters:
ji_st1 ((int, int)) – Index of the first InStamp.
ji_st2 ((int, int)) – Index of the second InStamp.
sim_mode (bool, optional) – Whether to count references without actually computing submatrices.
ji_st_out ((int, int)) – Index of the OutStamp. Needed for virtual memory.
- Returns:
The requested A submatrix. The shape is (st1.pix_cumsum[-1], st2.pix_cumsum[-1]).
- Return type:
np.array
- class SysMatB(blk)[source]
System matrix B attached to a coadd.Block instance.
The symtem matrix B is defined in Rowe+ 2011 Equation (18). IMPORTANT: The -2 coefficient is NOT included in this program.
- Parameters:
blk (pyimcom.coadd.Block) – The Block instance to which this SysMatB instance is attached.
- get_iosubmat(ji_st_in: int, int, ji_st_out: int, int, sim_mode: bool = False) numpy.array[source]
Return a requested B submatrix.
Note that this is the courterpart of the combination of SysMatA.get_iisubmat and the _compute_iisubmats method behind it, since SysMatB is much simpler than SysMatA.
- Parameters:
ji_st_in ((int, int)) – Index of the InStamp.
ji_st_out ((int, int)) – Index of the OutStamp.
sim_mode (bool, optional) – Whether to count references without actually computing submatrices.
- Returns:
np.array, shape – The requested B submatrix.
- Return type:
same as PSFOvl._call_io_cross
See also
SysMatA._compute_iisubmatsSee here for sim_mode description.