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.

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 Distributions
using Pigeons
using AstroImages
using CairoMakie

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

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

# Download sample images from GitHub
download(
    "https://github.com/sefffal/Octofitter.jl/raw/main/docs/image-examples-1.fits",
    "image-examples-1.fits"
)

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

5-element Vector{AstroImages.AstroImageMat{Float64, Tuple{DimensionalData.Dimensions.X{DimensionalData.Dimensions.Lookups.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, DimensionalData.Dimensions.Y{DimensionalData.Dimensions.Lookups.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}, Tuple{}, Matrix{Float64}, Tuple{DimensionalData.Dimensions.X{DimensionalData.Dimensions.Lookups.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}, DimensionalData.Dimensions.Y{DimensionalData.Dimensions.Lookups.Sampled{Int64, Base.OneTo{Int64}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Regular{Int64}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.NoMetadata}}}}}:
 [-0.340230578841865 -0.31239887526186355 … -0.2903753168231458 -0.34686827004814064; -0.29826568928566355 -0.27872689920999305 … -0.246452640103626 -0.2873140658333837; … ; 0.2673577868138141 0.2674003307584391 … -0.35696041486602914 -0.4320582965232348; 0.34147717271077094 0.3420908537177248 … -0.3358061981896538 -0.4139117427343716]
 [-1.0477769487233675 -0.903332487139985 … 1.4299602191511185 1.5070937596027922; -0.9888724404777535 -0.8589483767137767 … 1.1201264194935676 1.1803510698319304; … ; -0.43042419506343155 -0.46545013356425746 … -0.7066171674807686 -0.8619933383603304; -0.5837006039738425 -0.6021768042432977 … -0.8745923304009366 -1.0622653729194567]
 [-1.3673176596472436 -1.1985348376204175 … 1.1486692344915643 1.2929199629807906; -1.1143572848991925 -0.9810874969329584 … 0.9739105840286464 1.079311141942771; … ; 1.0292715351056276 0.8022139684685823 … 0.674123592556085 0.8519386886571362; 1.2170463982677187 0.9582307659215409 … 0.7602550494569255 0.959778160219427]
 [-0.16664741987337547 -0.13437529569419648 … -1.2963013699793633 -1.413565369547865; -0.2149990713151032 -0.1847368741065949 … -1.0549153738841368 -1.1490484959706464; … ; -0.5872462145987273 -0.5703179209210852 … 0.19469039731388912 0.31932314143450863; -0.7173113967617594 -0.6848676990867931 … 0.19122281859925802 0.3434506497145075]
 [0.47470453920043826 0.4747486393369815 … -0.16662161045865695 -0.16297895606178783; 0.4078627130832368 0.4041383409077344 … -0.18122203442704718 -0.17883512905931084; … ; 0.07003384675445126 0.19916060264117244 … -0.44137489323835993 -0.5253547777211444; 0.023766763613137065 0.1809119688664255 … -0.5699244785815155 -0.6776496918698217]

You can preview the image using imview from AstroImages:

# imshow2(image1, cmap=:magma) # for a single image
hcat(imview.(images, 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 observations:

image_dat = Table(;
    epoch = [1238.6, 1584.7, 3220.0, 7495.9, 7610.4],
    image = [
        AstroImages.recenter(images[1]),
        AstroImages.recenter(images[2]),
        AstroImages.recenter(images[3]),
        AstroImages.recenter(images[4]),
        AstroImages.recenter(images[5])
    ],
    platescale = [10.0, 10.0, 10.0, 10.0, 10.0]
)

image_like = ImageLikelihood(
    image_dat,
    name="SPHERE",
    variables=@variables begin
        # Planet flux in image units -- could be contrast, mags, Jy, or arb. as long as it's consistent
        flux ~ Normal(3.8, 0.5)
        # The following are optional parameters for marginalizing over instrument systematics:
        # Platescale uncertainty multiplier [could use: platescale ~ truncated(Normal(1, 0.01), lower=0)]
        platescale = 1.0
        # North angle offset in radians [could use: northangle ~ Normal(0, deg2rad(1))]
        northangle = 0.0
    end
)

OctofitterImages.ImageLikelihood Table with 5 columns and 5 rows:
     epoch   image                 platescale  contrast              ⋯
   ┌──────────────────────────────────────────────────────────────────
 1 │ 1238.6  [-0.340231 -0.31239…  10.0        65-element extrapol…  ⋯
 2 │ 1584.7  [-1.04778 -0.903332…  10.0        65-element extrapol…  ⋯
 3 │ 3220.0  [-1.36732 -1.19853 …  10.0        65-element extrapol…  ⋯
 4 │ 7495.9  [-0.166647 -0.13437…  10.0        65-element extrapol…  ⋯
 5 │ 7610.4  [0.474705 0.474749 …  10.0        65-element extrapol…  ⋯

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_b = Planet(
    name="b",
    basis=Visual{KepOrbit},
    likelihoods=[image_like],
    variables=@variables begin
        a ~ truncated(Normal(13, 4), lower=0.1, upper=100)
        e ~ Uniform(0.0, 0.5)
        i ~ Sine()
        M = system.M
        ω ~ UniformCircular()
        Ω ~ UniformCircular()
        θ ~ UniformCircular()
        tp = θ_at_epoch_to_tperi(θ, 1238.6; M, e, a, i, ω, Ω)
    end
)

Planet b
Derived:
  ω, Ω, θ, M, tp, 
Priors:
  a	Truncated(Distributions.Normal{Float64}(μ=13.0, σ=4.0); lower=0.1, upper=100.0)
  e	Distributions.Uniform{Float64}(a=0.0, b=0.5)
  i	Sine()
  ωx	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  ωy	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  Ωx	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  Ωy	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  θx	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  θy	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
OctofitterImages.ImageLikelihood Table with 5 columns and 5 rows:
     epoch   image                 platescale  contrast              ⋯
   ┌──────────────────────────────────────────────────────────────────
 1 │ 1238.6  [-0.340231 -0.31239…  10.0        65-element extrapol…  ⋯
 2 │ 1584.7  [-1.04778 -0.903332…  10.0        65-element extrapol…  ⋯
 3 │ 3220.0  [-1.36732 -1.19853 …  10.0        65-element extrapol…  ⋯
 4 │ 7495.9  [-0.166647 -0.13437…  10.0        65-element extrapol…  ⋯
 5 │ 7610.4  [0.474705 0.474749 …  10.0        65-element extrapol…  ⋯Octofitter.UnitLengthPrior{:ωx, :ωy}: √(ωx^2+ωy^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:Ωx, :Ωy}: √(Ωx^2+Ωy^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:θx, :θy}: √(θx^2+θy^2) ~ LogNormal(log(1), 0.02)

Note how we provided a prior on the photometry called flux in the variables block of the ImageLikelihood. This represents the expected flux of the planet in the image units.

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

Finally, create the system and pass in the planet.

sys = System(
    name="HD82134",
    companions=[planet_b],
    likelihoods=[],
    variables=@variables begin
        M ~ truncated(Normal(2.0, 0.1),lower=0.1)
        plx ~ truncated(Normal(45., 0.02),lower=0.1)
    end
)

System model HD82134
Derived:
  
Priors:
  M	Truncated(Distributions.Normal{Float64}(μ=2.0, σ=0.1); lower=0.1)
  plx	Truncated(Distributions.Normal{Float64}(μ=45.0, σ=0.02); lower=0.1)
Planet b
Derived:
  ω, Ω, θ, M, tp, 
Priors:
  a	Truncated(Distributions.Normal{Float64}(μ=13.0, σ=4.0); lower=0.1, upper=100.0)
  e	Distributions.Uniform{Float64}(a=0.0, b=0.5)
  i	Sine()
  ωx	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  ωy	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  Ωx	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  Ωy	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  θx	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
  θy	Distributions.Normal{Float64}(μ=0.0, σ=1.0)
OctofitterImages.ImageLikelihood Table with 5 columns and 5 rows:
     epoch   image                 platescale  contrast              ⋯
   ┌──────────────────────────────────────────────────────────────────
 1 │ 1238.6  [-0.340231 -0.31239…  10.0        65-element extrapol…  ⋯
 2 │ 1584.7  [-1.04778 -0.903332…  10.0        65-element extrapol…  ⋯
 3 │ 3220.0  [-1.36732 -1.19853 …  10.0        65-element extrapol…  ⋯
 4 │ 7495.9  [-0.166647 -0.13437…  10.0        65-element extrapol…  ⋯
 5 │ 7610.4  [0.474705 0.474749 …  10.0        65-element extrapol…  ⋯Octofitter.UnitLengthPrior{:ωx, :ωy}: √(ωx^2+ωy^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:Ωx, :Ωy}: √(Ωx^2+Ωy^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:θx, :θy}: √(θx^2+θy^2) ~ LogNormal(log(1), 0.02)

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

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(sys)

chain, pt = octofit_pigeons(model, n_rounds=10)
display(chain)

[ Info: Preparing model
┌ Info: Determined number of free variables
└   D = 12
┌ Info: Determined number type
└   T = Float64
ℓπcallback(θ): 0.000016 seconds
∇ℓπcallback(θ): 0.000028 seconds (1 allocation: 32 bytes)
┌ Warning: This model has priors that cannot be sampled IID.
└ @ OctofitterPigeonsExt ~/work/Octofitter.jl/Octofitter.jl/ext/OctofitterPigeonsExt/OctofitterPigeonsExt.jl:95
[ Info: Sampler running with multiple threads     : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: Performing global optimization with 12 parameters (0 initial parameter value provided).
┌ Warning: Unrecognized stop reason: Too many steps (101) without any function evaluations (probably search has converged). Defaulting to ReturnCode.Default.
└ @ Optimization ~/.julia/packages/Optimization/hfEs4/src/utils.jl:93
┌ Info: Found sample of initial positions
│   logpost_range = (22.99395417089064, 32.0413350848656)
└   mean_logpost = 28.69543252990779
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
  scans     restarts      Λ        Λ_var      time(s)    allc(B)  log(Z₁/Z₀)   min(α)     mean(α)    min(αₑ)   mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
        2          0       3.96       3.68      0.265   4.14e+06       13.8    0.00128      0.754          1          1
        4          0       3.44       4.65      0.118   1.39e+05       13.9   0.000475      0.739          1          1
        8          0       4.23       4.12      0.229    2.6e+05       12.8      0.351       0.73          1          1
       16          0       4.28       5.65      0.455   4.02e+05       11.3      0.044       0.68          1          1
       32          0       5.14       4.66      0.908   6.85e+05         14      0.284      0.684          1          1
       64          0          6       4.37       2.18   5.46e+07       24.9      0.133      0.665          1          1
      128          1        6.1       4.33          4   9.98e+07       55.3      0.357      0.664          1          1
      256         16        5.5       4.58       7.98   1.96e+08       64.9      0.317      0.675          1          1
      512         34       5.07       4.57       16.1   3.87e+08       64.9      0.309      0.689          1          1
 1.02e+03         59       5.28        4.6       31.7   7.75e+08       65.2      0.476      0.681          1          1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Diagnostics

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:

lines(
    chain["b_a"][:],
    axis=(;
        xlabel="iteration",
        ylabel="semi-major axis (aU)"
    )
)

And an auto-correlation plot:

using StatsBase
lines(
    autocor(chain["b_e"][:], 1:500),
    axis=(;
        xlabel="lag",
        ylabel="autocorrelation",
    )
)

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

Analysis

We can now view the orbit fit:

fig = octoplot(model, chain)

With a bit of work, we can plot one of the images under the orbit.

fig = octoplot(model, chain)
ax = fig.content[1] # grap first axis in the figure

# We have to do some annoying work to get the image orientated correctly,
# since we want the RA axis increasing to the left.
image_idx = 2
platescale = image_dat.platescale[image_idx]
img = AstroImages.recenter(AstroImage(collect(image_dat.image[image_idx])[end:-1:begin,:]))
imgax1 = dims(img,1) .* platescale
imgax2 = dims(img,2) .* platescale
h = heatmap!(ax, imgax1, imgax2, collect(img), colormap=:greys)
Makie.translate!(h, 0,0,-1) # Send heatmap to back of the plot

# Add colorbar for image
Colorbar(fig[1,2], h, label="image flux")

Makie.resize_to_layout!(fig)
fig

Another useful view would be the orbits over a stack of the maximum pixel values of all images.

fig = octoplot(model, chain)
ax = fig.content[1] # grap first axis in the figure

# We have to do some annoying work to get the image orientated correctly
# since we want the RA axis increasing to the left.
platescale = image_dat.platescale[image_idx]
imgs = maximum(stack(image_dat.image),dims=3)[:,:]
img = AstroImages.recenter(AstroImage(imgs[end:-1:begin,:]))
imgax1 = dims(img,1) .* platescale
imgax2 = dims(img,2) .* platescale
h = heatmap!(ax, imgax1, imgax2, collect(img), colormap=:greys)
Makie.translate!(h, 0,0,-1) # Send heatmap to back of the plot
Makie.resize_to_layout!(fig)
fig

Pair Plot

We can show the relationships between variables on a pair plot (aka corner plot):

using CairoMakie, PairPlots
octocorner(model, chain, small=true)

Note that this time, we also show the recovered photometry in the 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:

hist(chain["b_SPHERE_flux"][:], axis=(xlabel="flux", ylabel="counts"))

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

flux = chain["b_SPHERE_flux"]
snr = mean(flux)/std(flux)

14.711433937196686

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