bug fix regarding subpixel localization. Added stacking and drizzle

.gitignore

@@ -1,2 +1,4 @@
 Manifest.toml
 docs/build/
+examples/example_data/*
+examples/example_data/

Project.toml

@@ -13,13 +13,15 @@ ImageTransformations = "02fcd773-0e25-5acc-982a-7f6622650795"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 PSFModels = "9ba017d1-7760-46cd-84a3-1e79e9ae9ddc"
 Photometry = "af68cb61-81ac-52ed-8703-edc140936be4"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 Combinatorics = "1"
 CoordinateTransformations = "0.6"
 Distances = "0.10"
 ImageTransformations = "0.10"
 NearestNeighbors = "0.4"
-PSFModels = "0.8.2"
+PSFModels = "0.8"
 Photometry = "0.9"
+StaticArrays = "1.9"
 julia = "1.10"

examples/Project.toml

@@ -0,0 +1,11 @@
+[deps]
+AstroImages = "fe3fc30c-9b16-11e9-1c73-17dabf39f4ad"
+Astroalign = "7f4629bd-323a-4fad-ad42-4ee56350f27f"
+Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
+FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
+MultifileArrays = "18a1d154-b244-4181-98dd-8f3786468fa3"
+NDTools = "98581153-e998-4eef-8d0d-5ec2c052313d"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+View5D = "90d841e0-6953-4e90-9f3a-43681da8e949"

notebook.jl → examples/notebook.jl

@@ -20,7 +20,7 @@ end
 begin
 	import Pkg
 	Pkg.resolve()
-	Pkg.activate(joinpath(@__DIR__, "docs"))
+	Pkg.activate(joinpath(dirname(@__DIR__), "docs"))
 	Pkg.instantiate()
 
 	using Revise

examples/preprocess_helpers.jl

@@ -0,0 +1,22 @@
+using NDTools: select_region_view
+
+"""
+    correct_dark_flat(data, dark_img=nothing, flat_img=nothing, channel=nothing; T=Float32)
+
+corrects data using a `dark_img` and `flat_img` by subtracting the dark and deviding by the normalized, dark-subtracted `flat_img`.
+
+* minval: a minum value for the dark-subtracted flat image to avoid division by zero.
+"""
+function correct_dark_flat(data, dark_img=nothing, flat_img=nothing, channel=nothing, minval=0.001f0; T=Float32)
+    if !isnothing(dark_img)
+        dark_img = select_region_view(dark_img, new_size=size(data)[1:2]);
+        data = T.(data) .- T.(dark_img);
+    end
+    if !isnothing(flat_img)
+        dark_img = isnothing(dark_img) ? 0 : dark_img
+        flat_img = max.(select_region_view(flat_img, new_size=size(data)[1:2]) .- dark_img, minval);
+        flat_img ./= T.(sum(flat_img)/length(flat_img))        
+        data = T.(T.(data) ./ flat_img);
+    end
+    data
+end

examples/stack_measurement.jl

@@ -0,0 +1,71 @@
+#  This script stacks measured data and presents the result
+
+using Astroalign
+# using Pluto
+using Plots
+using Images # to display the stacked result
+# using View5D # for interactive viewing
+using Statistics: mean, median
+using NDTools: select_region
+using AstroImages # to be able to load FITS data
+using FileIO
+using MultifileArrays: load_series # add https://github.com/JuliaIO/MultifileArrays.jl.git
+
+include("preprocess_helpers.jl") # for background subtraction and field flatness correction
+
+if (false) # run this code to download the data and store it locally
+    # Here is some example data taken on my dwarf 3 telescope:
+    using Downloads
+    mydir = mktempdir()
+    zip_url = "https://cloud.uni-jena.de/s/9dqYsosX2DxTT5G/download"
+    # download the data (takes some time, about 500 Mb!)
+    Downloads.download(zip_url, joinpath(mydir, "data.zip"))
+    # unzip the downloaded data (in Windows only?). This creates the folder examples/example_data
+    run(`tar -xf $(joinpath(mydir, "data.zip"))`)  # works on Windows with tar
+end
+
+folder = "example_data\\"
+file_dark15 = raw"dark_exp_15.000000_gain_60_bin_1_12C_stack_9.fits"
+file_flat = raw"flat_gain_2_bin_1_ir_0.fits"
+files = raw"M 35_15s60_Astro_20260215-*_16C.fits"
+
+# data = load_series(load, files);
+dark = load(joinpath(folder, file_dark15))
+flat = load(joinpath(folder, file_flat))
+cd("example_data") # No idea why Windows cannot deal with the file path?
+data = load_series(load, files);
+cd("..")
+data = correct_dark_flat(data, dark, flat); # conversion to Float32 seems essential. In Float64 the fits seem to fail! 
+
+box_size = (15, 15)  
+ap_radius = 0.6 * first(box_size);
+stacked_d, all_params_d = stack_many_drizzle(data; dist_limit = 2, # use only one triangle
+        box_size, ap_radius, min_sigma = 1.5, nsigma = 1, ref_slice = 1, drizzle_supersampling = 2.0);
+
+prepare_for_viewer(v) = sqrt.(max.(0, reshape(v .- median(v, dims=(1,2)), (size(v)[1:2]...,1,size(v,3)))))
+prepare_for_display(v, m=10.0) = m .*colorview(RGB,permutedims(prepare_for_viewer(v), (4, 2, 1, 3))[:,:,:,1])
+# display the result as an RGB image
+
+plot(prepare_for_display(stacked_d, 0.08)) 
+# @vt prepare_for_viewer(stacked_d)
+
+# bin first and process the binned color data
+all_binned = Astroalign.bin_rgb(data);
+box_size = (9, 9)
+ap_radius = 0.6 * first(box_size);
+stacked_c, all_params_c = stack_many(all_binned; dist_limit = 2, 
+        box_size, ap_radius, min_sigma = 2.0, nsigma = 0.5, ref_slice = 1);
+
+# @vt prepare_for_viewer(stacked_c)
+plot(prepare_for_display(stacked_c, 0.08)) 
+
+# bin and sum colors first and process the binned monochrome data
+all_binned_m = Astroalign.bin_mono(data);
+box_size = (9, 9)
+ap_radius = 0.6 * first(box_size);
+stacked_m, all_param_m = stack_many(all_binned_m; dist_limit = 2, 
+        box_size, ap_radius, min_sigma = 2.0, nsigma = 0.5, ref_slice = 1);
+
+plot(prepare_for_display(cat(stacked_m,stacked_m,stacked_m, dims=3), 0.08)) 
+# heatmap(sqrt.(clamp.(stacked_m[:,:,1,1], 200, 250)))
+@vt prepare_for_viewer(stacked_m)

src/Astroalign.jl

@@ -3,191 +3,18 @@ module Astroalign
 using Combinatorics: combinations
 using CoordinateTransformations: kabsch
 using Distances: euclidean
-using ImageTransformations: warp
+using ImageTransformations: compose, warp, AffineMap, Translation
 using NearestNeighbors: nn, KDTree
 using PSFModels: gaussian, fit
-using Photometry: estimate_background,
-                  extract_sources,
-                  photometry,
-                  sigma_clip,
-                  CircularAperture,
-                  PeakMesh
+using Photometry: estimate_background, extract_sources, photometry, sigma_clip, CircularAperture, PeakMesh
+using StaticArrays: SVector
 
-export align_frame, get_sources, find_nearest, triangle_invariants
+export align_frame, find_nearest, get_sources, stack_many, stack_many_drizzle, triangle_invariants
 
-
-"""
-    get_sources(img; box_size = nothing, nsigma = 1)
-
-Extract candidate sources in `img` according to [`Photometry.Detection.extract_sources`](@extref). By default, `img` is first sigma clipped and then background subtracted before the candidate sources are extracted. `box_size` is passed to [`Photometry.Background.estimate_background`](@extref), and `nsigma` is passed to [`Photometry.Detection.extract_sources`](@extref). See the [Photometry.jl](@extref) documentation for more.
-
-TODO: Pass more options to clipping, background estimating, and extraction methods in [Photometry.jl](@extref).
-"""
-function get_sources(img; box_size = _compute_box_size(img), nsigma = 1, N_max = 10)
-    # Background subtract `img`
-    clipped = sigma_clip(img, 1, fill = NaN)
-    bkg, bkg_rms = estimate_background(clipped, box_size)
-    subt = img .- bkg[axes(img)...]
-
-    return (
-        # Sort detected sources from brightest to darkest
-        first(extract_sources(PeakMesh(; box_size, nsigma), subt, bkg, true), N_max),
-        # And also return the inputs, handy for debugging and data viz
-        subt,
-        bkg,
-    )
-end
-
-Base.@kwdef struct PSF
-    model = gaussian
-    params = (;)
-    func_kwargs = (;)
-    kwargs = (; x_abstol = 2e-6)
-end
-
-(p::PSF)(img) = fit_psf(img, p)
-
-function fit_psf(img_ap, p)
-    # Normalize
-    psf_data = collect(Float32, img_ap)
-    psf_data ./= maximum(psf_data)
-
-    # Set params
-    y, x = if !hasproperty(p.params, :x) && !hasproperty(p.params, :y)
-        Tuple(argmax(psf_data))
-    else
-        p.params.x, p.params.y
-    end
-
-    fwhm = if !hasproperty(p.params, :fwhm)
-        _compute_box_size(img_ap)
-    else
-        p.params.fwhm
-    end
-
-    # TODO: Generalize to other PSFs
-    params = (; x, y, fwhm)
-
-    # Fit
-    psf_params, psf_model = fit(p.model, params, psf_data; func_kwargs=p.func_kwargs, p.kwargs...)
-
-    return (; psf_params, psf_model, psf_data)
-end
-
-# Internal function used by `align`
-# Calls to `Photometry.photometry` with reasonable defaults
-function _photometry(img, box_size, ap_radius, min_fwhm, nsigma, f; filter_fwhm)
-    # Sources, background subtracted image, background
-    sources, subt, _ = get_sources(img; box_size, nsigma)
-
-    # Define apertures
-    aps = CircularAperture.(sources.y, sources.x, ap_radius)
-
-    phot = photometry(aps, subt; f)
-
-    if filter_fwhm
-        filter!(phot) do source
-            hypot(min_fwhm...) ≤ hypot(source.aperture_f.psf_params.fwhm...)
-        end
-    end
-
-    sort!(phot; by = x -> hypot(x.aperture_f.psf_params.fwhm...), rev = true)
-
-    return phot
-end
-
-"""
-    triangle_invariants(phot)
-
-Returns all combinations (``C``) of three candidate point sources from the table of sources `phot` returned by [`Photometry.Aperture.photometry`](@extref), and the computed invariant ``\\mathscr M`` for each according to Eq. 3 from [_Beroiz, M., Cabral, J. B., & Sanchez, B. (2020)_](https://ui.adsabs.harvard.edu/abs/2020A%26C....3200384B/abstract).
-"""
-function triangle_invariants(phot)
-    C = combinations(phot, 3)
-    ℳ = map(C) do (pa, pb, pc)
-        a, b, c = (
-            (pa.ycenter, pa.xcenter),
-            (pb.ycenter, pb.xcenter),
-            (pc.ycenter, pc.xcenter),
-        )
-        Ls = sort!([euclidean(a, b), euclidean(b, c), euclidean(a, c)])
-        (Ls[3] / Ls[2], Ls[2] / Ls[1])
-    end |> stack
-    return C, ℳ
-end
-
-"""
-    find_nearest(C_to, ℳ_to, C_from, ℳ_from)
-
-Return the closes pair of three points between the `from` and `to` frames in the invariant ``\\mathscr M`` space as computed by [`Astroalign.triangle_invariants`](@ref).
-"""
-function find_nearest(C_to, ℳ_to, C_from, ℳ_from)
-    idxs, dists = nn(KDTree(ℳ_to), ℳ_from)
-    idx_from = argmin(dists)
-    idx_to = idxs[idx_from]
-    sol_to = collect(C_to)[idx_to]
-    sol_from = collect(C_from)[idx_from]
-    return sol_to, sol_from
-end
-
-"""
-    function align_frame(img_to, img_from;
-        [box_size],
-        [ap_radius],
-        [f],
-        [min_fwhm],
-        [nsigma],
-    )
-
-Align `img_from` onto `img_to`. See below for keyword arguments currently available to control this process.
-
-* `box_size`: The size of the grid cells (in pixels) used to extract candidate point sources to use for alignment. Defaults to a tenth of the greatest common denominator of the dimensions of `img_to`. See [Photometry.jl > Source Detection Algorithms](@extref Photometry Source-Detection-Algorithms) for more.
-* `ap_radius`: The radius of the apertures (in pixel) to place around each point source. Defaults to 60% of `first(box_size)`. See [Photometry.jl > Aperture Photometry](@extref Photometry Aperture-Photometry) for more.
-* `f`: The function to compute within each aperture. Defaults to a 2D Gaussian fitted to the aperture center. See the [Source characterization](https://juliaastro.org/Astroalign.jl/notebook.html#Source-characterization) section of the accompanying Pluto.jl notebook for more.
-* `min_fwhm`: The minimum FWHM (in pixels) that an extracted point source must have to be considered as a control point. Defaults to a fifth of the width of the first image. See [PSFModels.jl > Fitting data](@extref PSFModels Fitting-data) for more.
-* `nsigma`: The number of standard deviations above the estimated background that a source must be to be considered as a control point. Defaults to 1. See [Photometry.jl > Source Detection Algorithms](@extref Photometry Source-Detection-Algorithms) for more.
-"""
-function align_frame(img_to, img_from;
-        box_size = _compute_box_size(img_to),
-        ap_radius = 0.6 * first(box_size),
-        f = PSF(),
-        min_fwhm = box_size .÷ 5,
-        nsigma = 1,
-    )
-    # Step 1: Identify control points
-    phot_to = _photometry(img_to, box_size, ap_radius, min_fwhm, nsigma, f; filter_fwhm=true)
-    phot_from = _photometry(img_from, box_size, ap_radius, min_fwhm, nsigma, f; filter_fwhm=true)
-
-    # Step 2: Calculate invariants
-    C_to, ℳ_to = triangle_invariants(phot_to)
-    C_from, ℳ_from = triangle_invariants(phot_from)
-
-    # Step 3: Select nearest
-    sol_to, sol_from = find_nearest(C_to, ℳ_to, C_from, ℳ_from)
-
-    # Transform
-    point_map = map(sol_to, sol_from) do source_to, source_from
-        [source_from.xcenter, source_from.ycenter] => [source_to.xcenter, source_to.ycenter]
-    end
-
-    tfm = kabsch(last.(point_map) => first.(point_map); scale=false)
-
-    return (
-        warp(img_from, tfm, axes(img_to)),
-        (;
-            point_map,
-            tfm,
-            C_to,
-            ℳ_to,
-            C_from,
-            ℳ_from,
-       )
-    )
-end
-
-function _compute_box_size(img)
-    w = gcd(size(img)...) ÷ 10
-    box_width = iseven(w) ? w + 1 : w
-    return (box_width, box_width)
-end
+include("utils.jl")
+include("findpeaks.jl")
+include("register.jl")
+include("warp.jl")
+include("stacker.jl")
 
 end # module