# PSFModels.jl

## Installation

PSFModels can be added from the Julia package manager

julia>]

(@v1.6) pkg> add PSFModels

## Getting Started

To import the library

julia> using PSFModels

None of the models are exported to avoid namespace clashes, but it can be verbose to continuously rewrite PSFModels. You can either import names directly

julia> using PSFModels: gaussian

julia> model = gaussian(x=0, y=0, fwhm=8)

or you can create an alias for PSFModels

# julia version 1.5 or below
using PSFModels
const M = PSFModels
# julia version 1.6 or above
import PSFModels as M

model = M.gaussian(x=0, y=0, fwhm=10)
PSFModels.PSFModelsModule

PSFModels

Statistical models for constructing point-spread functions (PSFs).

Models

The following models are currently implemented

Parameters

In general, the PSFs have a position, a full-width at half-maximum (FWHM) measure, and an amplitude. The position follows a 1-based pixel coordinate system, where (1, 1) represents the center of the bottom left pixel. This matches the indexing style of Julia as well as DS9, IRAF, SourceExtractor, and WCS. The FWHM is a consistent scale parameter for the models. That means a gaussian with a FWHM of 5 will be visually similar to an airydisk with a FWHM of 5. All models support a scalar (isotropic) FWHM and a FWHM for each axis (diagonal), as well as arbitrarily rotating the PSF.

Pixel-convention

The pixel convention adopted here is that the bottom-left pixel center is (1, 1). The column-major memory layout of julia puts the x axis as the rows of a matrix and the y axis as the columns. In other words, the axes unpack like

xs, ys = axes(image)

Usage

Evaluating models

Directly evaluating the functions is the most straightforward way to use this package

julia> gaussian(0, 0; x=0, y=0, fwhm=3)
1.0

julia> gaussian(BigFloat, 0, 0; x=0, y=0, fwhm=3, amp=0.1, bkg=1)
1.100000000000000088817841970012523233890533447265625

We also provide "curried" versions of the functions, which allow you to specify the parameters and evaluate the PSF later

julia> model = gaussian(x=0, y=0, fwhm=3);

julia> model(0, 0)
1.0

If we want to collect the model into a dense matrix, simply iterate over indices

julia> inds = CartesianIndices((-2:2, -2:2));

5×5 Matrix{Float64}:
0.0850494  0.214311  0.291632  0.214311  0.0850494
0.214311   0.54003   0.734867  0.54003   0.214311
0.291632   0.734867  1.0       0.734867  0.291632
0.214311   0.54003   0.734867  0.54003   0.214311
0.0850494  0.214311  0.291632  0.214311  0.0850494

This makes it very easy to evaluate the PSF on the same axes as an image (array)

julia> img = randn(5, 5);

julia> model.(CartesianIndices(img))
5×5 Matrix{Float64}:
0.54003      0.214311     0.0459292    0.00531559   0.000332224
0.214311     0.0850494    0.018227     0.00210949   0.000131843
0.0459292    0.018227     0.00390625   0.000452087  2.82555e-5
0.00531559   0.00210949   0.000452087  5.2322e-5    3.27013e-6
0.000332224  0.000131843  2.82555e-5   3.27013e-6   2.04383e-7

this is trivially expanded to fit "stamps" in images

julia> big_img = randn(1000, 1000);

julia> stamp_inds = (750:830, 400:485);

julia> stamp = @view big_img[stamp_inds...];

julia> stamp_model = model.(CartesianIndices(stamp_inds));

or we can create a loss function for fitting PSFs without allocating any memory. We are simply iterating over the image array!

julia> using Statistics

julia> mse = mean(I -> (big_img[I] - model(I))^2, CartesianIndices(stamp_inds));

Fitting data

There exists a simple, yet powerful, API for fitting data with PSFModels.fit.

# fit is not exported to avoid namespace clashes
using PSFModels: fit

stamp_inds = # optionally choose indices to "cutout"

# use an isotropic Gaussian
params, synthpsf = fit(gaussian, (x=12, y=13, fwhm=3.2, amp=0.1),
data, stamp_inds)
# elliptical, rotated Gaussian
params, synthpsf = fit(gaussian, (x=12, y=13, fwhm=(3.2, 3.2), amp=0.1, theta=0)
data, stamp_inds)
# obscured Airy disk
params, synthpsf = fit(airydisk, (x=12, y=13, fwhm=3.2, amp=0.1, ratio=0.3),
data, stamp_inds)
# bivariate Moffat with arbitrary alpha
params, synthpsf = fit(moffat, (x=12, y=13, fwhm=(3.2, 3.2), amp=0.1, alpha=1),
data, stamp_inds)

Plotting

Finally, we provide plotting recipes (psfplot/psfplot!) from RecipesBase.jl, which can be seen in use in the API/Reference section.

using Plots

model = gaussian(x=0, y=0, fwhm=(8, 10), theta=12)
psfplot(model, -30:30, -30:30, colorbar_scale=:log10)
source

## Contributing and Support

If you would like to contribute, feel free to open a pull request. If you want to discuss something before contributing, head over to discussions and join or open a new topic. If you're having problems with something, please open an issue.