Fast, allocation-free point-spread function (PSF) representations
gaussian(ornormal)airydiskmoffat
From the Julia REPL
julia> ]
(@v1.5) pkg> add PSFModelsTo import the library
julia> using PSFModelsjulia> 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(fwhm=10)For more in-depth usage and examples, please see the documentation.
First, load the package
julia> using PSFModelsDirectly 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)
0.1000000000000000055511151231257827021181583404541015625We 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.0If we want to collect the model into a dense matrix, simply iterate over indices
julia> inds = CartesianIndices((-2:2, -2:2));
julia> model.(inds) # broadcasting
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.0850494This 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-7this 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));There exists a simple, yet powerful, API for fitting data with these PSF models. See the full documentation for more details and examples.
# `fit` is not exported to avoid namespace clashes
using PSFModels: fit
data = # load data
stamp_inds = # optionally choose indices to "cutout"
# use an isotropic Gaussian
P0 = (x=12, y=13, fwhm=3.2, amp=0.1)
params, synthpsf = fit(gaussian, P0, data, stamp_inds)
# elliptical, rotated Gaussian
P0 = (x=12, y=13, fwhm=(3.2, 3.2), amp=0.1, theta=0)
params, synthpsf = fit(gaussian, P0, data, stamp_inds)
# obscured Airy disk
P0 = (x=12, y=13, fwhm=3.2, amp=0.1, ratio=0.3)
params, synthpsf = fit(airydisk, P0, data, stamp_inds)
# bivariate Moffat with arbitrary alpha
P0 = (x=12, y=13, fwhm=(3.2, 3.2), amp=0.1, alpha=1)
# fixed ("frozen") rotation angle
func_kwargs = (;theta=15)
params, synthpsf = fit(moffat, P0, data, stamp_inds; func_kwargs)We provide simple user recipes from RecipesBase.jl, which can be called with psfplot/psfplot!
using Plots
inds = (1:30, 1:30)
model = airydisk(x=12, y=13, fwhm=(4.5, 6.7), theta=12, ratio=0.3)
psfplot(model, inds, colorbar_scale=:log10)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.