Fitting Images

One of the key features of Octofitter.jl is the ability to search for planets directly from images of the system. Sampling from images is much more computationally demanding than sampling from astrometry, but it allows for a few very powerful results:

  1. You can search for a planet that is not well detected in a single image

By this, we mean you can feed in images of a system with no clear detections, and see if a planet is hiding in the noise based off of its Kepelerian motion.

  1. Not detecting a planet in a given image can be almost as useful as a detection for constraining its orbit.

If you have a clear detection in one epoch, but no detection in another, Octofitter can use the image from the second epoch to rule out large swathes of possible orbits.

Sampling from images can be freely combined with any known astrometry points, as well as astrometric acceleration. See advanced models for more details.


Image modelling is supported in Octofitter via the extension package OctofitterImages. To install it, run pkg> add

Preparing images

The first step will be to load your images. For this, we will use our AstroImages.jl package.

Start by loading your images:

using Octofitter
using OctofitterImages
using AstroImages

# Load individual iamges
# image1 = load("image1.fits")
# image2 = load("image2.fits")

# Or slices from a cube:
# cube = load("cube1.fits")
# image1 = cube[:,:,1] 

# Or multi-extension FITS (this example)
images = load("image-examples-1.fits",:)

You can preview the image using imshow2 from AstroImages:

# imshow2(image1, cmap=:magma) # for a single image
], cmap=:magma, clims=(-1.0, 4.0))


Your images should either be convolved with a gaussian of diameter one λ/D, or be matched filtered. This is so that the values of the pixels in the image represent the photometry at that location.

If you want to perform the convolution in Julia, see ImageFiltering.jl.

Build the model

First, we create a table of our image data that will be attached to the Planet:

image_data = Images(
    (band=:H, image=AstroImages.recenter(images[1]), platescale=10.0, epoch=1238.6),
    (band=:H, image=AstroImages.recenter(images[2]), platescale=10.0, epoch=1584.7),
    (band=:H, image=AstroImages.recenter(images[3]), platescale=10.0, epoch=3220.0),
    (band=:H, image=AstroImages.recenter(images[4]), platescale=10.0, epoch=7495.9),
    (band=:H, image=AstroImages.recenter(images[5]), platescale=10.0, epoch=7610.4),

Provide one entry for each image you want to sample from. Ensure that each image has been re-centered so that index [0,0] is the position of the star. Areas of the image where there is no data should be filled with NaN and will not contribute to the likelihood of your model. platescale should be the pixel scale of your images, in milliarseconds / pixel. epoch should be the Modified Julian Day (MJD) that your image was taken. You can use the mjd("2021-09-09") function to calculate this for you. band should be a symbol that matches the name you supplied when you created the Planet.

By default, the contrast of the images is calculated automatically, but you can supply your own contrast curve as well by also passing contrast=OctofitterImages.contrast_interp(AstroImages.recenter(my_image)).

You can freely mix and match images from different instruments as long as you specify the correct platescale. You can also provide images from multiple bands and they will be sampled independently. If you wish to tie them together, see Connecting Mass with Photometry.

Now specify the planet:

@planet X Visual{KepOrbit} begin
    a ~ Normal(13, 3)
    e ~ TruncatedNormal(0.2, 0.2, 0, 1.0)
    τ ~ Normal(0.5, 1)
    ω ~ Normal(0.1, deg2rad(30.))
    i ~ Normal(0.6, deg2rad(10.))
    Ω ~ Normal(0.0, deg2rad(30.))
    H ~ Normal(3.8, 0.5)
end image_data

Note how we also provided a prior on the photometry called H. We can put any name we want here, as long as it's used consistently throughout the model specification.

See Fit AstrometryLikelihood for a description of the different orbital parameters, and conventions used.

Finally, create the system and pass in the planet.

@system HD82134 begin
    M ~ Normal(2.0, 0.1),
    plx ~ Normal(45., 0.02),
end X

If you want to search for two or more planets in the same images, just create multiple Planets and pass the same images to each. You'll need to adjust the priors in some way to prevent overlap.

You can also do some very clever things like searching for planets that are co-planar and/or have a specific resonance between their periods. To do this, put the planet of the system or base period as variables of the system and derive the planet variables from those values of the system.


Sampling from images is much more challenging than relative astrometry or proper motion anomaly, so the fitting process tends to take longer.

This is because the posterior is much "bumpier" with images. One way this manifests is very high tree depths. You might see a sampling report that says max_tree_depth_frac = 0.9 or even 1.0. To encourage the sampler to take larger steps and explore the images, it's recommended to lower the target acceptance ratio to around 0.5±0.2 and also increase the number of adapataion steps.

model = Octofitter.LogDensityModel(HD82134; autodiff=:ForwardDiff, verbosity=4) # defaults are ForwardDiff, and verbosity=0

chain = octofit(
    model, 0.75,
    adaptation =   8_000,
    iterations =  10_000,
    max_depth =      14,

Sampling directly from images is somewhat slower than from astrometry. This example takes roughly 7 minutes on my laptop.


The first thing you should do with your results is check a few diagnostics to make sure the sampler converged as intended.

The acceptance rate should be somewhat lower than when fitting just astrometry, e.g. around the 0.6 target.

You can make a trace plot:

    ylabel="semi-major axis (aU)"

And an auto-correlation plot:

using StatsBase
    autocor(chain["X_e"], 1:500),

For this model, there is somewhat higher correlation between samples. Some thinning to remove this correlation is recommended. autocorrelation plot


You can plot the model as usual:

using Plots


In this case, the model is shown overtop a stack of the input images to help you visualize which peaks contributed to the fit. The images are stacked using the maximum function, so that bright spots from all images appear at once. The colour scale is inverted, so that the brightest peaks are shown in black.

You can also specify a lims=1000 parameter to set limits of the images to +/- 1000 mas, as in this example.

Pair Plot

We can show the relationships between variables on a pair plot (aka corner plot) using PairPlots.jl

using CairoMakie, PairPlots
table = (;
    a=         chain["X_a"],
    H=         chain["X_H"],
    e=         chain["X_e"],
    τ=         chain["X_τ"],
    :a => "a (au)",
    :H => "H (arb.)",
    :e => "e ",
    :i => "i (\\degree)",
    :Ω => "\\Omega (\\degree)",
    :ω => "\\omega (\\degree)",
    :τ => "\\tau ",
pairplot(table; labels)

Note that this time, we also show the recovered photometry in the corner plot.

corner plot

Assessing Detections

To assess a detection, we can treat all the orbital variables as nuisance parameters. We start by plotting the marginal distribution of the flux parameter, H:

histogram(chain["X_H"], xlabel="H", label="")

corner plot

We can calculate an analog of the traditional signal to noise ratio (SNR) using that same histogram:

flux = chain["X_H"]
snr = mean(flux)/std(flux) # 13.35 in this example

It might be better to consider a related measure, like the median flux over the interquartile distance. This will depend on your application.