GitHub

Aligning Astronomical Images

Credit: Beroiz, M., Cabral, J. B., & Sanchez, B. (2020)

Motivation

Aligning images comes up a lot in astronomy, like for co-adding exposures or timeseries photometry. The problem is that it can be computationally expensive to accomplish this via the traditional plate solving approach where we first need to calculate the WCS coordinates in each frame via a routine like astrometry.net, and then perform the relevant coordinate transformations from there.

Enter astroalign.py. This neat Python package sidesteps all of this by directly matching common star patterns between images to build a point-to-point correspondence. This page outlines the Julia reimplementation packaged as JuliaAstro/Astroalign.jl.

How It Works

Beroiz, Cabral, & Sanchez use the fact that triangles can be uniquely characterised to match sets of three stars (asterisms) between images. This point-to-point correspondence then gives everything needed to compute the affine transformation.

For this implementation they use the invariant $\mathscr M$ (the pair of two independent ratios of a triangle's side lengths, $L_i$) to define this unique characterization:

\[\begin{align} &\mathscr M = (L_3/L_2,\ L_2/L_1), \\ &\text{where}\ L_3 > L_2 > L_1\;. \end{align}\]

Astroalign.jl accomplishes this in the following steps:

  1. Identify the N_max brightest point-like sources in img_from and img_to.
  2. Calculate all triangular asterisms formed from these sources.
  3. Build a 2 × 3 × 2 × N array of candidate triangle-level correspondences by matching each from-triangle to its nearest to-triangle in the invariant $\mathscr M$ space. Vertices are assigned via a canonical ordering that is invariant under rotation, so the positional correspondence between matched triangles is geometrically consistent.
  4. Run RANSAC (Fischler & Bolles, 1981) on the triangle matches to robustly identify the largest set of mutually consistent correspondences ("inliers").
  5. Refine the transformation via the Kabsch/Umeyama least-squares algorithm applied to all vertex pairs from all inlier triangle matches.
  6. Warp img_from to the coordinates of img_to.

The rest of this page will show how to align simulated images using this technique.

Packages

Here is a summary of the packages we will use in this walkthrough:

# Main package and visualization packages
using Astroalign, CairoMakie
using AstroImages: AstroImage, Percent, Zscale, set_cmap!
set_cmap!(:cividis)

# Simulated star fields
using Random: Xoshiro
using LinearAlgebra: I, norm
using Rotations: RotMatrix2
using CoordinateTransformations: LinearMap, Translation
using ImageTransformations: Periodic, warp

Star field generator ✨

We'll start by creating a simulated star field to align on. For simplicity, we'll just create 12 Gaussian point sources placed randomly in a 300 x 300 grid with some noise over the whole image:

function gaussian(idx::CartesianIndex; x, y, fwhm, amp)
    σ = fwhm / (2 * sqrt(2 * log(2)))
    r, c = Tuple(idx)
    return amp * exp(-((r - x)^2 + (c - y)^2) / (2 * σ^2))
end

function generate_model(rng, model, params, inds)
    cartinds = CartesianIndices(inds)
    psf = model.(cartinds; params..., amp = 30_000)
    noise = rand(rng, 1000:3000, size(psf))
    return psf .+ noise
end
generate_model (generic function with 1 method)
img_to = let
    # Initial setup
    rng = Xoshiro(1)
    img_size = (1:300, 1:300)
    N_sources = 12
    FWHMs = [rand(rng, 1:10) for _ in 1:N_sources]
    pad = 10 # Minimum space between the stars (in pixels)
    positions_to = rand(rng, 1+pad:pad:300-pad, N_sources, 2)

    # Build stars
    map(zip(eachrow(positions_to), FWHMs)) do ((x, y), fwhm)
        generate_model(rng, gaussian, (; x, y, fwhm), img_size)
    end
# We wrap in an `AstroImage` to enable nice plotting recipes
end |> sum |> AstroImage
Example block output

We'll next create a transformed image that we would like to align back to the original:

# "Truth" values for transformation
const SCALE_0, ROT_0, TRANS_0 = 0.8, π/8, [10, 7]

img_from = let
    tfm = Translation(TRANS_0...) ∘
        LinearMap(RotMatrix2(ROT_0)) ∘
        LinearMap(SCALE_0 * I)

    warp(img_to, tfm, axes(img_to); fillvalue = Periodic())
end |> AstroImage
Example block output

In this particular case, img_from is i) scaled by a factor of 0.8, ii) rotated counter-clockwise by 22.5°, and iii) translated by [10, 7] pixels to arrive at img_to. We will now show how to align this image and recover these initial transformation parameters with Astroalign.jl.

First, here are some convenience functions that we will use to visualize our results:

# Global colorbar lims
const ZMIN, ZMAX = let
    lims = Percent(99.5).((img_to, img_from))
    minimum(first, lims), maximum(last, lims)
end

set_theme!(;
    Axis = (; aspect = DataAspect(), xticks = LinearTicks(4), yticks = LinearTicks(4)),
    Image = (; colorrange = (ZMIN, ZMAX), colormap = :cividis),
    # For default aperure plots
    Scatter = (;
        cycle = [], # Disable so that `color` is not overriden
        marker = Circle,
        markersize = 36, # Roughly the aperture size used
        markerspace = :data,
        color = :transparent,
        strokewidth = 2,
        strokecolor = :lightgreen,
    ),
)

function plot_pair(img_left, img_right;
    titles = ["img_left", "img_right"],
    colorrange = (ZMIN, ZMAX),
)
    fig = Figure(; size = (600, 300), figure_padding = (0, 0, 0, 0))

    ax_from, p_from = image(fig[1, 1], AstroImage(img_left); colorrange)
    ax_from.title = first(titles)

    ax_to, p_to = image(fig[1, 2], AstroImage(img_right); colorrange)
    ax_to.title = last(titles)
    hideydecorations!(ax_to)

    colsize!(fig.layout, 1, Aspect(1, 1.0))
    colsize!(fig.layout, 2, Aspect(1, 1.0))

    fig
end
plot_pair (generic function with 1 method)

Usage

We now use the exported align_frames function to align our image:

# Available options
opts_phot = (;
    box_size = (31, 31),
    ap_radius = 18.6,
    min_fwhm = (3, 3),
    nsigma = 1,
    f = Astroalign.com_psf,
    N_max = 10,
    use_fitpos = true,
);
opts_ransac = (; scale = true, ransac_threshold = 3.0);
opts_refinement = (; final_iters = 3, opts_ransac...)

# Align
arr_from_aligned = align_frames(img_from, img_to; opts_phot..., opts_refinement...);

# Visualize
plot_pair(arr_from_aligned, img_to; titles = ["img_from (aligned)", "img_to"])
Example block output

That's it! See the next section for a brief analysis on how well we did.

Recovered transformation

We can investigate the underlying transformation process by calling find_transform:

tfm_aligned, params_aligned = find_transform(img_from, img_to; opts_phot..., opts_refinement...);
(AffineMap([0.7390083156988436 -0.3064806662679341; 0.3064806662679341 0.7390083156988436], [10.073870117755433, 6.973263912109871]), (point_map = [[175.5252577413608, 257.45425110181094] => [60.90695903804799, 251.0202595497734], [90.52122374979801, 116.82381123908706] => [41.12110658240998, 120.97584088480511], [153.00109724559346, 104.66022532357749] => [91.02938765387015, 131.30134335378727], [281.54274953160694, 186.50772239069119] => [161.0031225640413, 231.0886499440785], [224.21096188871311, 47.826979492726764] => [161.07411034508587, 111.04524699796468], [126.89131684628056, 74.4587805873509] => [81.09355890608214, 100.87202231581347]], correspondences = [175.5252577413608 90.52122374979801 153.00109724559346; 257.45425110181094 116.82381123908706 104.66022532357749;;; 60.90695903804799 41.12110658240998 91.02938765387015; 251.0202595497734 120.97584088480511 131.30134335378727;;;; 175.5252577413608 153.00109724559346 281.54274953160694; 257.45425110181094 104.66022532357749 186.50772239069119;;; 60.90695903804799 91.02938765387015 161.0031225640413; 251.0202595497734 131.30134335378727 231.0886499440785;;;; 224.21096188871311 175.5252577413608 153.00109724559346; 47.826979492726764 257.45425110181094 104.66022532357749;;; 161.07411034508587 60.90695903804799 91.02938765387015; 111.04524699796468 251.0202595497734 131.30134335378727;;;; … ;;;; 126.89131684628056 122.65752549542285 99.05946551586688; 74.4587805873509 292.69640815823016 234.6834503362918;;; 81.09355890608214 2.3333925295358817 3.428398627952971; 100.87202231581347 260.944982642498 210.00415238097455;;;; 23.70064268503847 126.89131684628056 122.65752549542285; 279.58500057889296 74.4587805873509 292.69640815823016;;; 2.3333925295358817 41.12110658240998 60.90695903804799; 260.944982642498 120.97584088480511 251.0202595497734;;;; 122.65752549542285 23.70064268503847 99.05946551586688; 292.69640815823016 279.58500057889296 234.6834503362918;;; 91.02938765387015 41.12110658240998 81.09355890608214; 131.30134335378727 120.97584088480511 100.87202231581347], inlier_idxs = [1, 2, 3, 4, 8, 9, 10, 14, 15, 19, 29, 30, 31, 35, 36, 40, 50, 51, 55, 65], C_from = [(1, 3, 2), (1, 2, 4), (5, 1, 2), (1, 6, 2), (2, 1, 7), (8, 2, 1), (9, 2, 1), (3, 4, 1), (5, 1, 3), (6, 1, 3)  …  (5, 7, 6), (5, 8, 6), (5, 9, 6), (8, 5, 7), (5, 9, 7), (8, 5, 9), (8, 6, 7), (6, 9, 7), (8, 6, 9), (9, 8, 7)], ℳ_from = [1.0639669397439355 1.0135033481351452 … 1.0519553047537993 1.1379359796759485; 2.426368411983235 1.1945769440089773 … 2.1866855682091684 1.4006602356642557], C_to = [(3, 2, 1), (1, 4, 2), (5, 2, 1), (1, 2, 6), (7, 2, 1), (1, 8, 2), (1, 9, 2), (1, 10, 2), (4, 1, 3), (3, 1, 5)  …  (8, 6, 7), (6, 7, 9), (10, 6, 7), (8, 6, 9), (8, 6, 10), (6, 10, 9), (7, 8, 9), (8, 7, 10), (9, 7, 10), (10, 8, 9)], ℳ_to = [1.0482967802147551 1.5178096564155996 … 1.2727110267066146 1.1323061944899973; 1.6677355617223057 1.7000919004429984 … 1.4889115452483415 1.4835496715779886], phot_from = @NamedTuple{xcenter::Float64, ycenter::Float64, aperture_sum::Float64, aperture_f::@NamedTuple{psf_params::@NamedTuple{x::Float64, y::Float64, fwhm::Tuple{Float64, Float64}}, psf_model::String, psf_data::Photometry.Aperture.WeightedApertureCutout{Float64, Photometry.Aperture.CircularAperture{Float64}, AstroImageMat{Float64, Tuple{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}}, 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{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}}, 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{UnitRange{Int64}, UnitRange{Int64}}, typeof(*)}}}[(xcenter = 175.5252577413608, ycenter = 257.45425110181094, aperture_sum = 5.066811095980595e6, aperture_f = (psf_params = (x = 20.525257741360807, y = 18.45425110181096, fwhm = (9.598294444377542, 9.885632298094334)), psf_model = "com", psf_data = [0.0 -0.0 … -0.0 0.0; -0.0 -0.0 … 0.0 0.0; … ; -0.0 -0.0 … 0.0 0.0; -0.0 -0.0 … 0.0 0.0])), (xcenter = 153.00109724559346, ycenter = 104.66022532357749, aperture_sum = 3.3648372205258934e6, aperture_f = (psf_params = (x = 20.001097245593453, y = 19.66022532357749, fwhm = (8.187465091490306, 8.578817231091586)), psf_model = "com", psf_data = [0.0 0.0 … -0.0 -0.0; 0.0 0.0 … -0.0 -0.0; … ; -0.0 -0.0 … -0.0 -0.0; -0.0 -0.0 … -0.0 -0.0])), (xcenter = 90.52122374979801, ycenter = 116.82381123908706, aperture_sum = 3.4302831546280133e6, aperture_f = (psf_params = (x = 21.521223749798008, y = 20.823811239087068, fwhm = (8.09436635783948, 8.032014357017841)), psf_model = "com", psf_data = [-0.0 -0.0 … 0.0 -0.0; 0.0 0.0 … 0.0 0.0; … ; -0.0 -0.0 … 0.0 0.0; -0.0 -0.0 … 0.0 0.0])), (xcenter = 281.54274953160694, ycenter = 186.50772239069119, aperture_sum = 2.6173951166405436e6, aperture_f = (psf_params = (x = 19.542749531606965, y = 19.507722390691193, fwhm = (7.122246029681981, 7.228420617198278)), psf_model = "com", psf_data = [-0.0 -0.0 … -0.0 0.0; -0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … -0.0 -0.0; 0.0 0.0 … 0.0 -0.0])), (xcenter = 224.21096188871311, ycenter = 47.826979492726764, aperture_sum = 2.514023081816981e6, aperture_f = (psf_params = (x = 20.210961888713108, y = 19.82697949272676, fwhm = (7.105899759578943, 6.995090158867927)), psf_model = "com", psf_data = [-0.0 -0.0 … -0.0 0.0; -0.0 -0.0 … -0.0 -0.0; … ; 0.0 0.0 … 0.0 -0.0; 0.0 0.0 … 0.0 -0.0])), (xcenter = 126.89131684628056, ycenter = 74.4587805873509, aperture_sum = 1.8121249610262665e6, aperture_f = (psf_params = (x = 18.891316846280564, y = 20.458780587350898, fwhm = (6.3468448105288635, 6.133989513631395)), psf_model = "com", psf_data = [-0.0 -0.0 … -0.0 -0.0; 0.0 0.0 … -0.0 -0.0; … ; 0.0 0.0 … -0.0 -0.0; 0.0 0.0 … 0.0 0.0])), (xcenter = 99.05946551586688, ycenter = 234.6834503362918, aperture_sum = 156414.53547701085, aperture_f = (psf_params = (x = 20.05946551586688, y = 19.68345033629179, fwhm = (5.8669904996914095, 5.982076688845445)), psf_model = "com", psf_data = [-0.0 -0.0 … 0.0 0.0; 0.0 -0.0 … -0.0 -0.0; … ; 0.0 -0.0 … 0.0 0.0; 0.0 0.0 … -0.0 -0.0])), (xcenter = 23.70064268503847, ycenter = 279.58500057889296, aperture_sum = 1.1677841484072998e6, aperture_f = (psf_params = (x = 19.70064268503847, y = 19.58500057889298, fwhm = (5.266789407143618, 5.3283501363993855)), psf_model = "com", psf_data = [-0.0 -0.0 … -0.0 -0.0; -0.0 0.0 … -0.0 0.0; … ; 0.0 0.0 … 0.0 -0.0; -0.0 0.0 … -0.0 -0.0])), (xcenter = 122.65752549542285, ycenter = 292.69640815823016, aperture_sum = 865430.6965371289, aperture_f = (psf_params = (x = 19.65752549542284, y = 19.696408158230167, fwhm = (4.33147933382505, 4.996944768239902)), psf_model = "com", psf_data = [0.0 0.0 … -0.0 -0.0; 0.0 -0.0 … 0.7272302088289421 -0.0; … ; -0.0 -0.0 … -77.44939372979997 -0.0; -0.0 -0.0 … 0.0 -0.0]))], phot_to = @NamedTuple{xcenter::Float64, ycenter::Float64, aperture_sum::Float64, aperture_f::@NamedTuple{psf_params::@NamedTuple{x::Float64, y::Float64, fwhm::Tuple{Float64, Float64}}, psf_model::String, psf_data::Photometry.Aperture.WeightedApertureCutout{Float64, Photometry.Aperture.CircularAperture{Float64}, AstroImageMat{Float64, Tuple{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}}, 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{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}}, 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{UnitRange{Int64}, UnitRange{Int64}}, typeof(*)}}}[(xcenter = 3.428398627952971, ycenter = 210.00415238097455, aperture_sum = 197721.84992221327, aperture_f = (psf_params = (x = 12.428398627952971, y = 19.004152380974546, fwhm = (11.649406087952395, 14.210831637546242)), psf_model = "com", psf_data = [0.0 0.0 … -0.0 -0.0; -0.0 -0.0 … -0.0 -0.0; … ; -0.0 -0.0 … 0.0 0.0; -0.0 -0.0 … 0.0 0.0])), (xcenter = 60.90695903804799, ycenter = 251.0202595497734, aperture_sum = 3.235300005975329e6, aperture_f = (psf_params = (x = 19.906959038047983, y = 20.020259549773407, fwhm = (8.309517772983494, 8.361799058838034)), psf_model = "com", psf_data = [0.0 0.0 … 0.0 -0.0; -0.0 0.0 … 0.0 0.0; … ; -0.0 0.0 … 0.0 0.0; 0.0 -0.0 … 0.0 0.0])), (xcenter = 91.02938765387015, ycenter = 131.30134335378727, aperture_sum = 2.354980277806911e6, aperture_f = (psf_params = (x = 20.02938765387014, y = 19.301343353787278, fwhm = (7.515036765806923, 8.518149676427571)), psf_model = "com", psf_data = [0.0 -0.0 … 0.0 -0.0; 0.0 0.0 … -0.0 -0.0; … ; 0.0 -0.0 … -0.0 -0.0; -0.0 -0.0 … -0.0 -0.0])), (xcenter = 180.95446365945273, ycenter = 251.04941102421884, aperture_sum = 1.4110046147113673e6, aperture_f = (psf_params = (x = 19.95446365945273, y = 21.049411024218852, fwhm = (7.40650367626323, 7.562277886633494)), psf_model = "com", psf_data = [0.0 0.0 … 0.0 -0.0; 0.0 0.0 … 0.0 -0.0; … ; -0.0 -0.0 … 0.0 0.0; 0.0 0.0 … -0.0 0.0])), (xcenter = 41.12110658240998, ycenter = 120.97584088480511, aperture_sum = 2.1472886464067777e6, aperture_f = (psf_params = (x = 21.121106582409983, y = 20.975840884805116, fwhm = (7.217123334041633, 7.581501889101911)), psf_model = "com", psf_data = [0.0 0.0 … -0.0 -0.0; 0.0 -0.0 … 0.0 -0.0; … ; 0.0 0.0 … 0.0 -0.0; -0.0 0.0 … -0.0 0.0])), (xcenter = 2.3333925295358817, ycenter = 260.944982642498, aperture_sum = 563069.1469334947, aperture_f = (psf_params = (x = 11.333392529535882, y = 19.94498264249802, fwhm = (6.649343610432951, 7.658859768149461)), psf_model = "com", psf_data = [0.0 -0.0 … 0.0 -0.0; 0.0 0.0 … 0.0 0.0; … ; -0.0 0.0 … -0.0 -0.0; -0.0 -0.0 … 0.0 0.0])), (xcenter = 81.09355890608214, ycenter = 100.87202231581347, aperture_sum = 1.2605090520336095e6, aperture_f = (psf_params = (x = 19.09355890608214, y = 19.87202231581346, fwhm = (6.775271357141324, 7.513676051675258)), psf_model = "com", psf_data = [0.0 0.0 … -0.0 -0.0; 0.0 0.0 … -0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 -0.0 … -0.0 -0.0])), (xcenter = 241.26245071087774, ycenter = 220.79807733262214, aperture_sum = 770583.8214322324, aperture_f = (psf_params = (x = 20.262450710877744, y = 19.79807733262213, fwhm = (7.0434492712523875, 7.0519719359952315)), psf_model = "com", psf_data = [-0.0 0.0 … 0.0 0.0; -0.0 -0.0 … 0.0 -0.0; … ; -0.0 0.0 … -0.0 -0.0; 0.0 0.0 … -0.0 0.0])), (xcenter = 161.0031225640413, ycenter = 231.0886499440785, aperture_sum = 1.4279039899237866e6, aperture_f = (psf_params = (x = 20.003122564041302, y = 19.088649944078483, fwhm = (6.303858315548082, 6.760738540182835)), psf_model = "com", psf_data = [0.0 0.0 … 0.0 0.0; 0.0 -0.0 … 0.0 0.0; … ; 0.0 -0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])), (xcenter = 161.07411034508587, ycenter = 111.04524699796468, aperture_sum = 1.4074609870572896e6, aperture_f = (psf_params = (x = 20.074110345085877, y = 20.045246997964682, fwhm = (6.6815179956713004, 6.364941917737438)), psf_model = "com", psf_data = [0.0 -0.0 … 0.0 -0.0; -0.0 -0.0 … -0.0 0.0; … ; -0.0 0.0 … -0.0 -0.0; 0.0 -0.0 … -0.0 -0.0]))], phot_from_params = (sources = [(x = 187, y = 282, value = 31910.75498392784), (x = 48, y = 224, value = 31611.78284105163), (x = 116, y = 89, value = 31448.85977579732), (x = 74, y = 128, value = 31288.059206095317), (x = 259, y = 175, value = 31265.755156632447), (x = 105, y = 153, value = 30495.368152969477), (x = 280, y = 24, value = 29658.540322244273), (x = 235, y = 99, value = 29072.38502866167), (x = 293, y = 123, value = 28679.79152651803)], subt = [-909.2415731757756 -644.6286661319682 … -559.3054958877547 350.22187388814564; -215.67324086020744 -146.79245636952692 … -792.4315466917142 93.7175164961809; … ; -2013.5980328627593 700.9721667843514 … 881.9871524256923 541.7819468876405; -1406.2909580323649 -345.52173762462553 … -63.135625717401126 -566.4463179899794], bkg = [24259.08519263267 24259.08519263267 … 23928.04349016369 23928.04349016369; 24259.08519263267 24259.08519263267 … 23928.04349016369 23928.04349016369; … ; 24212.53091791701 24212.53091791701 … 24306.503027242627 24306.503027242627; 24212.53091791701 24212.53091791701 … 24306.503027242627 24306.503027242627], bkg_rms = [1051.2554984178919 1051.2554984178919 … 1018.5302778375594 1018.5302778375594; 1051.2554984178919 1051.2554984178919 … 1018.5302778375594 1018.5302778375594; … ; 1077.6616812103377 1077.6616812103377 … 1071.9551128818398 1071.9551128818398; 1077.6616812103377 1077.6616812103377 … 1071.9551128818398 1071.9551128818398], aps = Photometry.Aperture.CircularAperture{Float64}[CircularAperture(282.0, 187.0, r=18.6), CircularAperture(224.0, 48.0, r=18.6), CircularAperture(89.0, 116.0, r=18.6), CircularAperture(128.0, 74.0, r=18.6), CircularAperture(175.0, 259.0, r=18.6), CircularAperture(153.0, 105.0, r=18.6), CircularAperture(24.0, 280.0, r=18.6), CircularAperture(99.0, 235.0, r=18.6), CircularAperture(123.0, 293.0, r=18.6)]), phot_to_params = (sources = [(x = 232, y = 161, value = 32900.6254939205), (x = 211, y = 11, value = 32730.599989841456), (x = 101, y = 82, value = 32461.33699806252), (x = 251, y = 61, value = 32258.940285642646), (x = 111, y = 161, value = 32144.994738026708), (x = 120, y = 40, value = 32123.72803973405), (x = 132, y = 91, value = 31958.694224143386), (x = 221, y = 241, value = 31322.157358645025), (x = 261, y = 11, value = 31129.62592109082), (x = 250, y = 181, value = 28627.510831911794)], subt = [-2417.1268844221086 -0.1268844221085601 … -784.2363977485948 881.7636022514052; 4694.873115577891 -3110.1268844221086 … 2162.7636022514052 -50.236397748594754; … ; -1305.2397660818751 718.7602339181249 … -437.7011331444737 271.2988668555263; -1473.2397660818751 -640.2397660818751 … 488.2988668555263 -3136.7011331444737], bkg = [23966.12688442211 23966.12688442211 … 24115.236397748595 24115.236397748595; 23966.12688442211 23966.12688442211 … 24115.236397748595 24115.236397748595; … ; 23939.239766081875 23939.239766081875 … 23677.701133144474 23677.701133144474; 23939.239766081875 23939.239766081875 … 23677.701133144474 23677.701133144474], bkg_rms = [1324.7810024521868 1324.7810024521868 … 1344.461636998835 1344.461636998835; 1324.7810024521868 1324.7810024521868 … 1344.461636998835 1344.461636998835; … ; 1286.5948794420888 1286.5948794420888 … 1360.9184886654202 1360.9184886654202; 1286.5948794420888 1286.5948794420888 … 1360.9184886654202 1360.9184886654202], aps = Photometry.Aperture.CircularAperture{Float64}[CircularAperture(161.0, 232.0, r=18.6), CircularAperture(11.0, 211.0, r=18.6), CircularAperture(82.0, 101.0, r=18.6), CircularAperture(61.0, 251.0, r=18.6), CircularAperture(161.0, 111.0, r=18.6), CircularAperture(40.0, 120.0, r=18.6), CircularAperture(91.0, 132.0, r=18.6), CircularAperture(241.0, 221.0, r=18.6), CircularAperture(11.0, 261.0, r=18.6), CircularAperture(181.0, 250.0, r=18.6)])))
Tip

This can be used in conjunction with apply_transform to apply the transformation found. This is what align_frames calls under the hood.

Decomposing it into scale (S), rotation (R), and translation (T) components then gives:

function decompose_tfm(tfm)
    M = tfm.linear
    S = sqrt(M'M)
    R = M * inv(S)
    T = tfm.translation
    return (; S, R, T)
end

# To compare with "truth" values
p_diff(x, x0) = round(100 * (x - x0) / x0; digits = 3)
p_diff(x::AbstractVector, x0::AbstractVector) = round(100 * norm(x - x0) / norm(x0); digits = 3)
print_diff(name, x, x0) = println(
    "$(name) : $(round.(x; digits = 4)) (truth: $(x0), error: $(p_diff(x, x0))%)"
)
print_diff (generic function with 1 method)
S, R, T = decompose_tfm(tfm_aligned)

params_tfm = (scale = S[1], rot = atan(R[2, 1], R[1, 1]), trans = T)

print_diff("Scale", params_tfm.scale, SCALE_0)
print_diff("Rotation", rad2deg(params_tfm.rot), rad2deg(ROT_0))
print_diff("Translation", params_tfm.trans, TRANS_0)
Scale : 0.8 (truth: 0.8, error: 0.005%)
Rotation : 22.5247 (truth: 22.5, error: 0.11%)
Translation : [10.0739, 6.9733] (truth: [10, 7], error: 0.644%)

Taking a look at our RANSAC pass, these final transformation values were determined from 10 out of 35 detected correspondences (28.6 %).

n_inliers = length(params_aligned.inlier_idxs)
n_total   = size(params_aligned.correspondences, 4)
println("RANSAC inliers: $n_inliers / $n_total ($(round(100*n_inliers/n_total; digits = 1))%)")
RANSAC inliers: 20 / 84 (23.8%)

The rest of this document will walk through how this is accomplished behind the scenes, and the different options that we can pass to align_frames.

Step 1: Identify control points

This step is done solely on the Photometry.jl side for both our img_from and img_to images, which Astroalign.jl calls with some reasonable defaults via Astroalign._photometry.

phot_to, phot_to_params = Astroalign._photometry(img_to; opts_phot...)
phot_from, phot_from_params = Astroalign._photometry(img_from; opts_phot...)
(@NamedTuple{xcenter::Float64, ycenter::Float64, aperture_sum::Float64, aperture_f::@NamedTuple{psf_params::@NamedTuple{x::Float64, y::Float64, fwhm::Tuple{Float64, Float64}}, psf_model::String, psf_data::Photometry.Aperture.WeightedApertureCutout{Float64, Photometry.Aperture.CircularAperture{Float64}, AstroImageMat{Float64, Tuple{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}}, 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{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}}, 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{UnitRange{Int64}, UnitRange{Int64}}, typeof(*)}}}[(xcenter = 175.5252577413608, ycenter = 257.45425110181094, aperture_sum = 5.066811095980595e6, aperture_f = (psf_params = (x = 20.525257741360807, y = 18.45425110181096, fwhm = (9.598294444377542, 9.885632298094334)), psf_model = "com", psf_data = [0.0 -0.0 … -0.0 0.0; -0.0 -0.0 … 0.0 0.0; … ; -0.0 -0.0 … 0.0 0.0; -0.0 -0.0 … 0.0 0.0])), (xcenter = 153.00109724559346, ycenter = 104.66022532357749, aperture_sum = 3.3648372205258934e6, aperture_f = (psf_params = (x = 20.001097245593453, y = 19.66022532357749, fwhm = (8.187465091490306, 8.578817231091586)), psf_model = "com", psf_data = [0.0 0.0 … -0.0 -0.0; 0.0 0.0 … -0.0 -0.0; … ; -0.0 -0.0 … -0.0 -0.0; -0.0 -0.0 … -0.0 -0.0])), (xcenter = 90.52122374979801, ycenter = 116.82381123908706, aperture_sum = 3.4302831546280133e6, aperture_f = (psf_params = (x = 21.521223749798008, y = 20.823811239087068, fwhm = (8.09436635783948, 8.032014357017841)), psf_model = "com", psf_data = [-0.0 -0.0 … 0.0 -0.0; 0.0 0.0 … 0.0 0.0; … ; -0.0 -0.0 … 0.0 0.0; -0.0 -0.0 … 0.0 0.0])), (xcenter = 281.54274953160694, ycenter = 186.50772239069119, aperture_sum = 2.6173951166405436e6, aperture_f = (psf_params = (x = 19.542749531606965, y = 19.507722390691193, fwhm = (7.122246029681981, 7.228420617198278)), psf_model = "com", psf_data = [-0.0 -0.0 … -0.0 0.0; -0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … -0.0 -0.0; 0.0 0.0 … 0.0 -0.0])), (xcenter = 224.21096188871311, ycenter = 47.826979492726764, aperture_sum = 2.514023081816981e6, aperture_f = (psf_params = (x = 20.210961888713108, y = 19.82697949272676, fwhm = (7.105899759578943, 6.995090158867927)), psf_model = "com", psf_data = [-0.0 -0.0 … -0.0 0.0; -0.0 -0.0 … -0.0 -0.0; … ; 0.0 0.0 … 0.0 -0.0; 0.0 0.0 … 0.0 -0.0])), (xcenter = 126.89131684628056, ycenter = 74.4587805873509, aperture_sum = 1.8121249610262665e6, aperture_f = (psf_params = (x = 18.891316846280564, y = 20.458780587350898, fwhm = (6.3468448105288635, 6.133989513631395)), psf_model = "com", psf_data = [-0.0 -0.0 … -0.0 -0.0; 0.0 0.0 … -0.0 -0.0; … ; 0.0 0.0 … -0.0 -0.0; 0.0 0.0 … 0.0 0.0])), (xcenter = 99.05946551586688, ycenter = 234.6834503362918, aperture_sum = 156414.53547701085, aperture_f = (psf_params = (x = 20.05946551586688, y = 19.68345033629179, fwhm = (5.8669904996914095, 5.982076688845445)), psf_model = "com", psf_data = [-0.0 -0.0 … 0.0 0.0; 0.0 -0.0 … -0.0 -0.0; … ; 0.0 -0.0 … 0.0 0.0; 0.0 0.0 … -0.0 -0.0])), (xcenter = 23.70064268503847, ycenter = 279.58500057889296, aperture_sum = 1.1677841484072998e6, aperture_f = (psf_params = (x = 19.70064268503847, y = 19.58500057889298, fwhm = (5.266789407143618, 5.3283501363993855)), psf_model = "com", psf_data = [-0.0 -0.0 … -0.0 -0.0; -0.0 0.0 … -0.0 0.0; … ; 0.0 0.0 … 0.0 -0.0; -0.0 0.0 … -0.0 -0.0])), (xcenter = 122.65752549542285, ycenter = 292.69640815823016, aperture_sum = 865430.6965371289, aperture_f = (psf_params = (x = 19.65752549542284, y = 19.696408158230167, fwhm = (4.33147933382505, 4.996944768239902)), psf_model = "com", psf_data = [0.0 0.0 … -0.0 -0.0; 0.0 -0.0 … 0.7272302088289421 -0.0; … ; -0.0 -0.0 … -77.44939372979997 -0.0; -0.0 -0.0 … 0.0 -0.0]))], (sources = [(x = 187, y = 282, value = 31910.75498392784), (x = 48, y = 224, value = 31611.78284105163), (x = 116, y = 89, value = 31448.85977579732), (x = 74, y = 128, value = 31288.059206095317), (x = 259, y = 175, value = 31265.755156632447), (x = 105, y = 153, value = 30495.368152969477), (x = 280, y = 24, value = 29658.540322244273), (x = 235, y = 99, value = 29072.38502866167), (x = 293, y = 123, value = 28679.79152651803)], subt = [-909.2415731757756 -644.6286661319682 … -559.3054958877547 350.22187388814564; -215.67324086020744 -146.79245636952692 … -792.4315466917142 93.7175164961809; … ; -2013.5980328627593 700.9721667843514 … 881.9871524256923 541.7819468876405; -1406.2909580323649 -345.52173762462553 … -63.135625717401126 -566.4463179899794], bkg = [24259.08519263267 24259.08519263267 … 23928.04349016369 23928.04349016369; 24259.08519263267 24259.08519263267 … 23928.04349016369 23928.04349016369; … ; 24212.53091791701 24212.53091791701 … 24306.503027242627 24306.503027242627; 24212.53091791701 24212.53091791701 … 24306.503027242627 24306.503027242627], bkg_rms = [1051.2554984178919 1051.2554984178919 … 1018.5302778375594 1018.5302778375594; 1051.2554984178919 1051.2554984178919 … 1018.5302778375594 1018.5302778375594; … ; 1077.6616812103377 1077.6616812103377 … 1071.9551128818398 1071.9551128818398; 1077.6616812103377 1077.6616812103377 … 1071.9551128818398 1071.9551128818398], aps = Photometry.Aperture.CircularAperture{Float64}[CircularAperture(282.0, 187.0, r=18.6), CircularAperture(224.0, 48.0, r=18.6), CircularAperture(89.0, 116.0, r=18.6), CircularAperture(128.0, 74.0, r=18.6), CircularAperture(175.0, 259.0, r=18.6), CircularAperture(153.0, 105.0, r=18.6), CircularAperture(24.0, 280.0, r=18.6), CircularAperture(99.0, 235.0, r=18.6), CircularAperture(123.0, 293.0, r=18.6)]))

This performs source extraction and source characterization of our images, storing the results in the phot_from_params and phot_to_params named tuples above. Here is a preview of the detected soures in each image:

fig = plot_pair(img_from, img_to; titles = ["img_from", "img_to"])

# Show apertures
scatter!(fig.content[1], phot_from_params.sources.y, phot_from_params.sources.x)
scatter!(fig.content[2], phot_to_params.sources.y, phot_to_params.sources.x)

fig
Example block output

Source extraction

astroalign.py uses sep under the hood for its source extraction. We use Photometry.Detection.extract_sources to pull out the regions around the brightest pixels. A FWHM estimate is then computed for each detected source (see the next section) so the brightest, well-resolved sources can be selected as control points and noise-like detections (hot pixels, artifacts, etc.) can be filtered out via min_fwhm. Source extraction is executed by Astroalign._get_sources.

In the future, additional photometry options may be added.

Source characterization

Photometry.Aperture.photometry automatically computes aperture sums and returns them in a nice table for us. We also pass a function, Astroalign.com_psf, to compute PSF statistics for each source from a fast, non-iterative center-of-mass estimate, and store them in the table as well.

In addition to the usual photometry fields returned, the aperture_f field contains a named tuple of PSF information computed by default with Astroalign.com_psf:

  • psf_params: Named tuple of x and y center (relative to the aperture) and per-axis fwhm of the source.
  • psf_model: A tag indicating which characterization method was used ("com" for center-of-mass).
  • psf_data: The underlying array of the data within the aperture being characterized.

A different characterization function can be supplied via the f keyword of find_transform — any callable that takes an aperture cutout and returns a named tuple with the fields above will work.

Below is a quick visual check of an observed point source and its center-of-mass fit:

function inspect_psf(phot; kwargs...)
    psf_data = phot.aperture_f.psf_data
    (; x, y, fwhm) = phot.aperture_f.psf_params

    println((; phot.xcenter, phot.ycenter))
    println(phot.aperture_f.psf_params)

    σx, σy = fwhm ./ (2 * sqrt(2 * log(2)))
    model = [exp(-((r - x)^2 / (2σx^2) + (c - y)^2 / (2σy^2)))
             for r in axes(psf_data, 1), c in axes(psf_data, 2)]

    obs = AstroImage(psf_data)
    psf = AstroImage(model)

    plot_pair(obs, psf; kwargs...)
end

inspect_psf(first(phot_from); colorrange = (0, 1), titles = ["Data", "Model"])
Example block output
Note

PSF centers are relative to the aperture, while xcenter and ycenter are relative to the whole image. Astroalign.jl performs the necessary conversions from the former to the latter in Astroalign.to_subpixel before reporting the final fitted values.

With our sources identified, we now turn to the next step in the alignment algorithm.

Step 2: Calculate invariants

This is done internally in align_frames, but the invariants $\mathscr M_i$ can also be exposed with Astroalign._triangle_invariants.

C_to, ℳ_to = Astroalign._triangle_invariants(phot_to)

# This can also be accessed through the named tuple
# returned by `find_transform`.
(; C_from, ℳ_from) = params_aligned
(point_map = [[175.5252577413608, 257.45425110181094] => [60.90695903804799, 251.0202595497734], [90.52122374979801, 116.82381123908706] => [41.12110658240998, 120.97584088480511], [153.00109724559346, 104.66022532357749] => [91.02938765387015, 131.30134335378727], [281.54274953160694, 186.50772239069119] => [161.0031225640413, 231.0886499440785], [224.21096188871311, 47.826979492726764] => [161.07411034508587, 111.04524699796468], [126.89131684628056, 74.4587805873509] => [81.09355890608214, 100.87202231581347]], correspondences = [175.5252577413608 90.52122374979801 153.00109724559346; 257.45425110181094 116.82381123908706 104.66022532357749;;; 60.90695903804799 41.12110658240998 91.02938765387015; 251.0202595497734 120.97584088480511 131.30134335378727;;;; 175.5252577413608 153.00109724559346 281.54274953160694; 257.45425110181094 104.66022532357749 186.50772239069119;;; 60.90695903804799 91.02938765387015 161.0031225640413; 251.0202595497734 131.30134335378727 231.0886499440785;;;; 224.21096188871311 175.5252577413608 153.00109724559346; 47.826979492726764 257.45425110181094 104.66022532357749;;; 161.07411034508587 60.90695903804799 91.02938765387015; 111.04524699796468 251.0202595497734 131.30134335378727;;;; … ;;;; 126.89131684628056 122.65752549542285 99.05946551586688; 74.4587805873509 292.69640815823016 234.6834503362918;;; 81.09355890608214 2.3333925295358817 3.428398627952971; 100.87202231581347 260.944982642498 210.00415238097455;;;; 23.70064268503847 126.89131684628056 122.65752549542285; 279.58500057889296 74.4587805873509 292.69640815823016;;; 2.3333925295358817 41.12110658240998 60.90695903804799; 260.944982642498 120.97584088480511 251.0202595497734;;;; 122.65752549542285 23.70064268503847 99.05946551586688; 292.69640815823016 279.58500057889296 234.6834503362918;;; 91.02938765387015 41.12110658240998 81.09355890608214; 131.30134335378727 120.97584088480511 100.87202231581347], inlier_idxs = [1, 2, 3, 4, 8, 9, 10, 14, 15, 19, 29, 30, 31, 35, 36, 40, 50, 51, 55, 65], C_from = [(1, 3, 2), (1, 2, 4), (5, 1, 2), (1, 6, 2), (2, 1, 7), (8, 2, 1), (9, 2, 1), (3, 4, 1), (5, 1, 3), (6, 1, 3)  …  (5, 7, 6), (5, 8, 6), (5, 9, 6), (8, 5, 7), (5, 9, 7), (8, 5, 9), (8, 6, 7), (6, 9, 7), (8, 6, 9), (9, 8, 7)], ℳ_from = [1.0639669397439355 1.0135033481351452 … 1.0519553047537993 1.1379359796759485; 2.426368411983235 1.1945769440089773 … 2.1866855682091684 1.4006602356642557], C_to = [(3, 2, 1), (1, 4, 2), (5, 2, 1), (1, 2, 6), (7, 2, 1), (1, 8, 2), (1, 9, 2), (1, 10, 2), (4, 1, 3), (3, 1, 5)  …  (8, 6, 7), (6, 7, 9), (10, 6, 7), (8, 6, 9), (8, 6, 10), (6, 10, 9), (7, 8, 9), (8, 7, 10), (9, 7, 10), (10, 8, 9)], ℳ_to = [1.0482967802147551 1.5178096564155996 … 1.2727110267066146 1.1323061944899973; 1.6677355617223057 1.7000919004429984 … 1.4889115452483415 1.4835496715779886], phot_from = @NamedTuple{xcenter::Float64, ycenter::Float64, aperture_sum::Float64, aperture_f::@NamedTuple{psf_params::@NamedTuple{x::Float64, y::Float64, fwhm::Tuple{Float64, Float64}}, psf_model::String, psf_data::Photometry.Aperture.WeightedApertureCutout{Float64, Photometry.Aperture.CircularAperture{Float64}, AstroImageMat{Float64, Tuple{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}}, 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{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}}, 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{UnitRange{Int64}, UnitRange{Int64}}, typeof(*)}}}[(xcenter = 175.5252577413608, ycenter = 257.45425110181094, aperture_sum = 5.066811095980595e6, aperture_f = (psf_params = (x = 20.525257741360807, y = 18.45425110181096, fwhm = (9.598294444377542, 9.885632298094334)), psf_model = "com", psf_data = [0.0 -0.0 … -0.0 0.0; -0.0 -0.0 … 0.0 0.0; … ; -0.0 -0.0 … 0.0 0.0; -0.0 -0.0 … 0.0 0.0])), (xcenter = 153.00109724559346, ycenter = 104.66022532357749, aperture_sum = 3.3648372205258934e6, aperture_f = (psf_params = (x = 20.001097245593453, y = 19.66022532357749, fwhm = (8.187465091490306, 8.578817231091586)), psf_model = "com", psf_data = [0.0 0.0 … -0.0 -0.0; 0.0 0.0 … -0.0 -0.0; … ; -0.0 -0.0 … -0.0 -0.0; -0.0 -0.0 … -0.0 -0.0])), (xcenter = 90.52122374979801, ycenter = 116.82381123908706, aperture_sum = 3.4302831546280133e6, aperture_f = (psf_params = (x = 21.521223749798008, y = 20.823811239087068, fwhm = (8.09436635783948, 8.032014357017841)), psf_model = "com", psf_data = [-0.0 -0.0 … 0.0 -0.0; 0.0 0.0 … 0.0 0.0; … ; -0.0 -0.0 … 0.0 0.0; -0.0 -0.0 … 0.0 0.0])), (xcenter = 281.54274953160694, ycenter = 186.50772239069119, aperture_sum = 2.6173951166405436e6, aperture_f = (psf_params = (x = 19.542749531606965, y = 19.507722390691193, fwhm = (7.122246029681981, 7.228420617198278)), psf_model = "com", psf_data = [-0.0 -0.0 … -0.0 0.0; -0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … -0.0 -0.0; 0.0 0.0 … 0.0 -0.0])), (xcenter = 224.21096188871311, ycenter = 47.826979492726764, aperture_sum = 2.514023081816981e6, aperture_f = (psf_params = (x = 20.210961888713108, y = 19.82697949272676, fwhm = (7.105899759578943, 6.995090158867927)), psf_model = "com", psf_data = [-0.0 -0.0 … -0.0 0.0; -0.0 -0.0 … -0.0 -0.0; … ; 0.0 0.0 … 0.0 -0.0; 0.0 0.0 … 0.0 -0.0])), (xcenter = 126.89131684628056, ycenter = 74.4587805873509, aperture_sum = 1.8121249610262665e6, aperture_f = (psf_params = (x = 18.891316846280564, y = 20.458780587350898, fwhm = (6.3468448105288635, 6.133989513631395)), psf_model = "com", psf_data = [-0.0 -0.0 … -0.0 -0.0; 0.0 0.0 … -0.0 -0.0; … ; 0.0 0.0 … -0.0 -0.0; 0.0 0.0 … 0.0 0.0])), (xcenter = 99.05946551586688, ycenter = 234.6834503362918, aperture_sum = 156414.53547701085, aperture_f = (psf_params = (x = 20.05946551586688, y = 19.68345033629179, fwhm = (5.8669904996914095, 5.982076688845445)), psf_model = "com", psf_data = [-0.0 -0.0 … 0.0 0.0; 0.0 -0.0 … -0.0 -0.0; … ; 0.0 -0.0 … 0.0 0.0; 0.0 0.0 … -0.0 -0.0])), (xcenter = 23.70064268503847, ycenter = 279.58500057889296, aperture_sum = 1.1677841484072998e6, aperture_f = (psf_params = (x = 19.70064268503847, y = 19.58500057889298, fwhm = (5.266789407143618, 5.3283501363993855)), psf_model = "com", psf_data = [-0.0 -0.0 … -0.0 -0.0; -0.0 0.0 … -0.0 0.0; … ; 0.0 0.0 … 0.0 -0.0; -0.0 0.0 … -0.0 -0.0])), (xcenter = 122.65752549542285, ycenter = 292.69640815823016, aperture_sum = 865430.6965371289, aperture_f = (psf_params = (x = 19.65752549542284, y = 19.696408158230167, fwhm = (4.33147933382505, 4.996944768239902)), psf_model = "com", psf_data = [0.0 0.0 … -0.0 -0.0; 0.0 -0.0 … 0.7272302088289421 -0.0; … ; -0.0 -0.0 … -77.44939372979997 -0.0; -0.0 -0.0 … 0.0 -0.0]))], phot_to = @NamedTuple{xcenter::Float64, ycenter::Float64, aperture_sum::Float64, aperture_f::@NamedTuple{psf_params::@NamedTuple{x::Float64, y::Float64, fwhm::Tuple{Float64, Float64}}, psf_model::String, psf_data::Photometry.Aperture.WeightedApertureCutout{Float64, Photometry.Aperture.CircularAperture{Float64}, AstroImageMat{Float64, Tuple{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}}, 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{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}}, 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{UnitRange{Int64}, UnitRange{Int64}}, typeof(*)}}}[(xcenter = 3.428398627952971, ycenter = 210.00415238097455, aperture_sum = 197721.84992221327, aperture_f = (psf_params = (x = 12.428398627952971, y = 19.004152380974546, fwhm = (11.649406087952395, 14.210831637546242)), psf_model = "com", psf_data = [0.0 0.0 … -0.0 -0.0; -0.0 -0.0 … -0.0 -0.0; … ; -0.0 -0.0 … 0.0 0.0; -0.0 -0.0 … 0.0 0.0])), (xcenter = 60.90695903804799, ycenter = 251.0202595497734, aperture_sum = 3.235300005975329e6, aperture_f = (psf_params = (x = 19.906959038047983, y = 20.020259549773407, fwhm = (8.309517772983494, 8.361799058838034)), psf_model = "com", psf_data = [0.0 0.0 … 0.0 -0.0; -0.0 0.0 … 0.0 0.0; … ; -0.0 0.0 … 0.0 0.0; 0.0 -0.0 … 0.0 0.0])), (xcenter = 91.02938765387015, ycenter = 131.30134335378727, aperture_sum = 2.354980277806911e6, aperture_f = (psf_params = (x = 20.02938765387014, y = 19.301343353787278, fwhm = (7.515036765806923, 8.518149676427571)), psf_model = "com", psf_data = [0.0 -0.0 … 0.0 -0.0; 0.0 0.0 … -0.0 -0.0; … ; 0.0 -0.0 … -0.0 -0.0; -0.0 -0.0 … -0.0 -0.0])), (xcenter = 180.95446365945273, ycenter = 251.04941102421884, aperture_sum = 1.4110046147113673e6, aperture_f = (psf_params = (x = 19.95446365945273, y = 21.049411024218852, fwhm = (7.40650367626323, 7.562277886633494)), psf_model = "com", psf_data = [0.0 0.0 … 0.0 -0.0; 0.0 0.0 … 0.0 -0.0; … ; -0.0 -0.0 … 0.0 0.0; 0.0 0.0 … -0.0 0.0])), (xcenter = 41.12110658240998, ycenter = 120.97584088480511, aperture_sum = 2.1472886464067777e6, aperture_f = (psf_params = (x = 21.121106582409983, y = 20.975840884805116, fwhm = (7.217123334041633, 7.581501889101911)), psf_model = "com", psf_data = [0.0 0.0 … -0.0 -0.0; 0.0 -0.0 … 0.0 -0.0; … ; 0.0 0.0 … 0.0 -0.0; -0.0 0.0 … -0.0 0.0])), (xcenter = 2.3333925295358817, ycenter = 260.944982642498, aperture_sum = 563069.1469334947, aperture_f = (psf_params = (x = 11.333392529535882, y = 19.94498264249802, fwhm = (6.649343610432951, 7.658859768149461)), psf_model = "com", psf_data = [0.0 -0.0 … 0.0 -0.0; 0.0 0.0 … 0.0 0.0; … ; -0.0 0.0 … -0.0 -0.0; -0.0 -0.0 … 0.0 0.0])), (xcenter = 81.09355890608214, ycenter = 100.87202231581347, aperture_sum = 1.2605090520336095e6, aperture_f = (psf_params = (x = 19.09355890608214, y = 19.87202231581346, fwhm = (6.775271357141324, 7.513676051675258)), psf_model = "com", psf_data = [0.0 0.0 … -0.0 -0.0; 0.0 0.0 … -0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 -0.0 … -0.0 -0.0])), (xcenter = 241.26245071087774, ycenter = 220.79807733262214, aperture_sum = 770583.8214322324, aperture_f = (psf_params = (x = 20.262450710877744, y = 19.79807733262213, fwhm = (7.0434492712523875, 7.0519719359952315)), psf_model = "com", psf_data = [-0.0 0.0 … 0.0 0.0; -0.0 -0.0 … 0.0 -0.0; … ; -0.0 0.0 … -0.0 -0.0; 0.0 0.0 … -0.0 0.0])), (xcenter = 161.0031225640413, ycenter = 231.0886499440785, aperture_sum = 1.4279039899237866e6, aperture_f = (psf_params = (x = 20.003122564041302, y = 19.088649944078483, fwhm = (6.303858315548082, 6.760738540182835)), psf_model = "com", psf_data = [0.0 0.0 … 0.0 0.0; 0.0 -0.0 … 0.0 0.0; … ; 0.0 -0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0])), (xcenter = 161.07411034508587, ycenter = 111.04524699796468, aperture_sum = 1.4074609870572896e6, aperture_f = (psf_params = (x = 20.074110345085877, y = 20.045246997964682, fwhm = (6.6815179956713004, 6.364941917737438)), psf_model = "com", psf_data = [0.0 -0.0 … 0.0 -0.0; -0.0 -0.0 … -0.0 0.0; … ; -0.0 0.0 … -0.0 -0.0; 0.0 -0.0 … -0.0 -0.0]))], phot_from_params = (sources = [(x = 187, y = 282, value = 31910.75498392784), (x = 48, y = 224, value = 31611.78284105163), (x = 116, y = 89, value = 31448.85977579732), (x = 74, y = 128, value = 31288.059206095317), (x = 259, y = 175, value = 31265.755156632447), (x = 105, y = 153, value = 30495.368152969477), (x = 280, y = 24, value = 29658.540322244273), (x = 235, y = 99, value = 29072.38502866167), (x = 293, y = 123, value = 28679.79152651803)], subt = [-909.2415731757756 -644.6286661319682 … -559.3054958877547 350.22187388814564; -215.67324086020744 -146.79245636952692 … -792.4315466917142 93.7175164961809; … ; -2013.5980328627593 700.9721667843514 … 881.9871524256923 541.7819468876405; -1406.2909580323649 -345.52173762462553 … -63.135625717401126 -566.4463179899794], bkg = [24259.08519263267 24259.08519263267 … 23928.04349016369 23928.04349016369; 24259.08519263267 24259.08519263267 … 23928.04349016369 23928.04349016369; … ; 24212.53091791701 24212.53091791701 … 24306.503027242627 24306.503027242627; 24212.53091791701 24212.53091791701 … 24306.503027242627 24306.503027242627], bkg_rms = [1051.2554984178919 1051.2554984178919 … 1018.5302778375594 1018.5302778375594; 1051.2554984178919 1051.2554984178919 … 1018.5302778375594 1018.5302778375594; … ; 1077.6616812103377 1077.6616812103377 … 1071.9551128818398 1071.9551128818398; 1077.6616812103377 1077.6616812103377 … 1071.9551128818398 1071.9551128818398], aps = Photometry.Aperture.CircularAperture{Float64}[CircularAperture(282.0, 187.0, r=18.6), CircularAperture(224.0, 48.0, r=18.6), CircularAperture(89.0, 116.0, r=18.6), CircularAperture(128.0, 74.0, r=18.6), CircularAperture(175.0, 259.0, r=18.6), CircularAperture(153.0, 105.0, r=18.6), CircularAperture(24.0, 280.0, r=18.6), CircularAperture(99.0, 235.0, r=18.6), CircularAperture(123.0, 293.0, r=18.6)]), phot_to_params = (sources = [(x = 232, y = 161, value = 32900.6254939205), (x = 211, y = 11, value = 32730.599989841456), (x = 101, y = 82, value = 32461.33699806252), (x = 251, y = 61, value = 32258.940285642646), (x = 111, y = 161, value = 32144.994738026708), (x = 120, y = 40, value = 32123.72803973405), (x = 132, y = 91, value = 31958.694224143386), (x = 221, y = 241, value = 31322.157358645025), (x = 261, y = 11, value = 31129.62592109082), (x = 250, y = 181, value = 28627.510831911794)], subt = [-2417.1268844221086 -0.1268844221085601 … -784.2363977485948 881.7636022514052; 4694.873115577891 -3110.1268844221086 … 2162.7636022514052 -50.236397748594754; … ; -1305.2397660818751 718.7602339181249 … -437.7011331444737 271.2988668555263; -1473.2397660818751 -640.2397660818751 … 488.2988668555263 -3136.7011331444737], bkg = [23966.12688442211 23966.12688442211 … 24115.236397748595 24115.236397748595; 23966.12688442211 23966.12688442211 … 24115.236397748595 24115.236397748595; … ; 23939.239766081875 23939.239766081875 … 23677.701133144474 23677.701133144474; 23939.239766081875 23939.239766081875 … 23677.701133144474 23677.701133144474], bkg_rms = [1324.7810024521868 1324.7810024521868 … 1344.461636998835 1344.461636998835; 1324.7810024521868 1324.7810024521868 … 1344.461636998835 1344.461636998835; … ; 1286.5948794420888 1286.5948794420888 … 1360.9184886654202 1360.9184886654202; 1286.5948794420888 1286.5948794420888 … 1360.9184886654202 1360.9184886654202], aps = Photometry.Aperture.CircularAperture{Float64}[CircularAperture(161.0, 232.0, r=18.6), CircularAperture(11.0, 211.0, r=18.6), CircularAperture(82.0, 101.0, r=18.6), CircularAperture(61.0, 251.0, r=18.6), CircularAperture(161.0, 111.0, r=18.6), CircularAperture(40.0, 120.0, r=18.6), CircularAperture(91.0, 132.0, r=18.6), CircularAperture(241.0, 221.0, r=18.6), CircularAperture(11.0, 261.0, r=18.6), CircularAperture(181.0, 250.0, r=18.6)]))

Below is a plot comparing the compents of the computed invariants for all control points in our from and to images. Overlapping regions between the from and to clouds indicate similar triangles found by Astroalign.jl. Compare to Fig. 1 in Beroiz et al. (2020).

# Use default theme
with_theme() do
    fig, ax, p = scatter(ℳ_to[1, :], ℳ_to[2, :];
        label = "img_to",
    )

    ax.xlabel = "L3/L2"
    ax.ylabel = "L2/L1"

    scatter!(ax, ℳ_from[1, :], ℳ_from[2, :];
        markersize = 20,
        color = :transparent,
        strokewidth = 2,
        strokecolor = :cornflowerblue,
        label = "img_from",
    )

    axislegend()

    fig
end
Example block output
Note

The number of triangle combinations may differ between frames if sources drift towards or off the edge of the frame between images. All that is needed is one matching triangle.

Step 3: Build candidate correspondences

We next build our list of candidate correspondences in this invariant space via a nearest neighbors search:

correspondences = Astroalign._build_correspondences(C_from, ℳ_from, phot_from, C_to, ℳ_to, phot_to)

println("Candidate triangle matches: $(size(correspondences, 4))")
Candidate triangle matches: 84

Step 4: RANSAC pass

The largest mutually consistent set of correspondences ("inliers") is found via a RANSAC pass using JuliaAstro/ConsensusFitting.jl with Astroalign._ransac:

fwd_tfm_initial, inlier_idxs_initial = Astroalign._ransac(correspondences; opts_ransac...)

println("Initial RANSAC inliers: $(length(inlier_idxs_initial)) / $(size(correspondences, 4))")
Initial RANSAC inliers: 20 / 84

Step 5: Refine transformation

The transformation and inlier set from the previous step are successively refined via Astroalign._refine_transform using all detected control points, capturing previously missed inliers while dropping incorrect assignments:

tfm, inlier_idxs, point_map = Astroalign._refine_transform(fwd_tfm_initial, inlier_idxs_initial, correspondences; opts_refinement...)
(AffineMap([0.7390083156988436 -0.3064806662679341; 0.3064806662679341 0.7390083156988436], [10.073870117755433, 6.973263912109871]), [1, 2, 3, 4, 8, 9, 10, 14, 15, 19, 29, 30, 31, 35, 36, 40, 50, 51, 55, 65], [[175.5252577413608, 257.45425110181094] => [60.90695903804799, 251.0202595497734], [90.52122374979801, 116.82381123908706] => [41.12110658240998, 120.97584088480511], [153.00109724559346, 104.66022532357749] => [91.02938765387015, 131.30134335378727], [281.54274953160694, 186.50772239069119] => [161.0031225640413, 231.0886499440785], [224.21096188871311, 47.826979492726764] => [161.07411034508587, 111.04524699796468], [126.89131684628056, 74.4587805873509] => [81.09355890608214, 100.87202231581347]])

For this example, all 10 initial inliers remain after refinement:

println("Final RANSAC inliers: $(length(inlier_idxs)) / $(size(correspondences, 4))")

inlier_idxs == inlier_idxs_initial

The matched control points in both images are shown below:

fig = plot_pair(img_from, img_to; titles = ["img_from", "img_to"])

# Solution apertures
x_from, y_from, x_to, y_to = (getindex.(f.(point_map), i) for f in (first, last) for i in (1, 2))
strokecolor = Makie.categorical_colors(:tab10, length(x_from))
scatter!(fig.content[1], x_from, y_from; strokecolor)
scatter!(fig.content[2], x_to, y_to; strokecolor)

fig
Example block output

Step 6: Apply transformation

Once the linear transformation parameters have been finalized in Step 5, we hand it off to ImageTransformations.warp to perform resampling to align img_from with img_to:

img_aligned_from = AstroImage(warp(img_from, inv(tfm), axes(img_to)))

plot_pair(img_aligned_from, img_to; titles = ["img_from (aligned)", "img_to"])
Example block output