GitHub

Examples

Fitting a PSF

Here is a brief example which shows how to construct a loss function for fitting a PSFModel to some data.

using PSFModels
using PSFModels: fit
using HCIDatasets: BetaPictoris
using Plots
using Statistics

# convenience function for plotting
function imshow(data; kwargs...)
    xlim = extrema(axes(data, 1))
    ylim = extrema(axes(data, 2))
    heatmap(transpose(data); xlim=xlim, ylim=ylim,
            aspect_ratio=1, clims=(1e-5, Inf), kwargs...)
end

# get a PSF from HCIDatasets.jl;
# you may be prompted to download the file
psf = BetaPictoris[:psf]
inds = CartesianIndices(psf)

imshow(psf)
Example block output

We can fit this data with a variety of models, here showcasing the flexible PSFModels.fit function.

Gaussian

Using gaussian

params = (x=20, y=20, fwhm=5, amp=0.1)
P_gauss, mod_gauss = fit(gaussian, params, psf)
pairs(P_gauss)
pairs(::NamedTuple) with 4 entries:
  :x    => 20.0191
  :y    => 20.01
  :fwhm => 4.80293
  :amp  => 0.100216
plot(
    imshow(psf, title="Data"),
    imshow(mod_gauss.(inds), title="Model"),
    cbar=false,
    ticks=false,
    layout=2,
    size=(600, 300)
)
Example block output

and now using a rotated, elliptical Gaussian

params = (x=20, y=20, fwhm=(5, 5), amp=0.1, theta=0)
P_ellip, mod_ellip = fit(gaussian, params, psf)
pairs(P_ellip)
pairs(::NamedTuple) with 5 entries:
  :x     => 20.0183
  :y     => 20.01
  :fwhm  => (4.76706, 4.83663)
  :amp   => 0.100243
  :theta => 0.0
plot(
    imshow(psf, title="Data"),
    imshow(mod_ellip.(inds), title="Model"),
    cbar=false,
    ticks=false,
    layout=2,
    size=(600, 300)
)
Example block output

Airy disk

Now with airydisk

params = (x=20, y=20, fwhm=5, amp=0.1, ratio=0.3)
P_airy, mod_airy = fit(airydisk, params, psf)
pairs(P_airy)
pairs(::NamedTuple) with 5 entries:
  :x     => 19.9891
  :y     => 20.0062
  :fwhm  => 5.05739
  :amp   => 0.139425
  :ratio => 0.375189
plot(
    imshow(psf, title="Data"),
    imshow(mod_airy.(inds), title="Model"),
    cbar=false,
    ticks=false,
    layout=2,
    size=(600, 300)
)
Example block output

Moffat

And finally, with moffat

params = (x=20, y=20, fwhm=(5, 5), amp=0.1, theta=0, alpha=2)
P_moff, mod_moff = fit(moffat, params, psf)
pairs(P_moff)
pairs(::NamedTuple) with 6 entries:
  :x     => 20.0366
  :y     => 20.0154
  :fwhm  => (3.95567, 3.78329)
  :amp   => 0.11638
  :theta => 42.3939
  :alpha => 1.46268
plot(
    imshow(psf, title="Data"),
    imshow(mod_moff.(inds), title="Model"),
    cbar=false,
    ticks=false,
    layout=2,
    size=(600, 300)
)
Example block output

Changing optimization parameters

Any keyword arguments get passed on to Optim.optimize, and you can change the algorithm used with the alg keyword

# load Optim.jl to use the Newton method
using Optim

params = (x=20, y=20, fwhm=(5, 5), amp=0.1, theta=0, alpha=2)
P_moff, mod_moff = fit(moffat, params, psf; alg=Newton())
pairs(P_moff)
pairs(::NamedTuple) with 6 entries:
  :x     => 20.0359
  :y     => 20.0151
  :fwhm  => (3.87008, 3.8665)
  :amp   => 0.116352
  :theta => 0.0
  :alpha => 1.46154

We can also "freeze" parameters by creating a named tuple and passing it to func_kwargs

params = (;x=10, y=20, fwhm=(5, 5), amp=0.1)
func_kwargs = (;alpha=2)
P_moff2, mod_moff2 = fit(moffat, params, psf; func_kwargs)
pairs(P_moff2)
pairs(::NamedTuple) with 4 entries:
  :x    => 20.0349
  :y    => 20.0144
  :fwhm => (4.19137, 4.19478)
  :amp  => 0.110515
plot(
    imshow(psf, title="Data"),
    imshow(mod_moff2.(inds), title="Model"),
    cbar=false,
    ticks=false,
    layout=2,
    size=(600, 300)
)
Example block output