commit 92ed176ba250755bdfec88536921f1d0f69f0113 from: ale date: Mon Aug 17 04:41:06 2026 UTC Improve two intentional-design limitations: relax estimate_psf's isolation filter before falling back; warn when zogy_subtract's V_ast is silently skipped - estimate_psf: on an empty stamp list, halve min_separation and retry (relaxation_attempts, default 2) before falling back to the analytic Moffat fit — a real empirical PSF from fewer, closer stars still beats a parametric approximation, which the previous hard binary threw away. - zogy_subtract: warn (once per session) when both n_sources/r_sources are left unset, so omitting V_ast is a visible choice for a direct caller instead of a silent default they could miss. run_pipeline always supplies both, so this never fires on the normal pipeline path. Both are real, tested improvements to items previously documented as "intentional design, not a bug" — the design itself was still the right call in each case, but neither was as good as it could be within that design. Full narrative in docs/src/investigation-log.md, including a real fit_moffat_psf stamp-overlap limitation surfaced while writing the regression tests. commit - 8ff5933c85c68fcb18d5be2ee05c05d4d2c6df00 commit + 92ed176ba250755bdfec88536921f1d0f69f0113 blob - 413f5933626675d76786f401e56e3c9361681dbf blob + dd66e303cdccc4f3354ebc934d1667928a28cfc4 --- docs/src/index.md +++ docs/src/index.md @@ -108,18 +108,26 @@ all fixed with regression tests, and how each was diag ## Known limitations -- **Empirical PSF's quality depends on the field.** `estimate_psf` stacks - real star cutouts, which captures the true PSF shape (wings included) - without fitting a model family per instrument, but needs enough bright, - isolated, unsaturated stars to do it. When a field doesn't have them, it - falls back (by default) to `fit_moffat_psf` — a parametric Moffat fit — - rather than failing outright; the fallback trades exact PSF shape for - robustness, and is not itself a substitute for a genuinely well-behaved - field. -- **`zogy_subtract`'s astrometric-noise term (`V_ast`) is opt-in at the - `zogy_subtract` level** — it needs `n_sources`/`r_sources` passed - explicitly, and is `0` without them. `run_pipeline` always supplies - them, so this only matters when calling `zogy_subtract` directly. +- **Empirical PSF's quality still depends on the field, though less than + it used to.** `estimate_psf` stacks real star cutouts, which captures + the true PSF shape (wings included) without fitting a model family per + instrument, but needs enough bright, isolated, unsaturated stars to do + it. It now retries at a progressively relaxed `min_separation` (halved, + up to `relaxation_attempts` times, default 2) before giving up — a + field with a few usable stars at a tighter isolation radius still gives + the real PSF shape, which the analytic fallback never can. Only once + even the most relaxed attempt finds nothing does it fall back (by + default) to `fit_moffat_psf` — a parametric Moffat fit — rather than + failing outright; the fallback trades exact PSF shape for robustness, + and is not itself a substitute for a genuinely well-behaved field. +- **`zogy_subtract`'s astrometric-noise term (`V_ast`) is still opt-in at + the `zogy_subtract` level** — it needs `n_sources`/`r_sources` passed + explicitly, and is `0` without them; this is deliberate API layering + (`run_pipeline` always supplies them, so a direct caller who doesn't + need the extra `detect_sources` cost can skip it), not something + planned to change. What did change: leaving both unset now emits a + `@warn` (once per session), so omitting `V_ast` is a visible choice + instead of a silent default a direct caller could miss. - **`find_variable_sources` still has a real, measured false-positive floor on real single-epoch aperture photometry, even after fixing it once.** The first hypothesis tried — pixel-grid jitter in blob - 276447b0b4cc24d5b2fe594094385a85dc9f2bc8 blob + dc62fc44f1207d65635045e46ac2e212d7be5467 --- docs/src/investigation-log.md +++ docs/src/investigation-log.md @@ -468,3 +468,50 @@ real, independently-confirmed variable star invisible false-positive rate, while leaving the real variable's signal a 3x margin over threshold — the largest floor checked that doesn't cost real detections, not the smallest false-positive rate achievable. + +## Revisiting the two intentional-design "limitations" + +Two items in `docs/src/index.md`'s Known Limitations were design +decisions, not bugs — but "intentional" isn't the same as "as good as it +can be." Asked directly whether either could be genuinely improved +without abandoning the design choice behind it. + +**`estimate_psf`'s empirical-vs-fallback split was a hard binary — a +field either had stars passing the isolation/saturation filter at exactly +the given `min_separation`, or it fell all the way back to an analytic +Moffat fit.** But a moderately (not severely) crowded field might have +real, usable stars at a *slightly* tighter isolation radius — the current +code was throwing that away and jumping straight to an approximation +instead of trying harder for the real thing. Added `relaxation_attempts` +(default 2): on an empty stamp list, halve `min_separation` and retry, +up to that many times, before falling back. A real empirical PSF from +fewer, closer stars still beats a parametric approximation, which is the +whole reason `estimate_psf` exists over just always using +`fit_moffat_psf`. Regression-tested with two synthetic Gaussian sources +25 px apart (fails the default `min_separation=40`, but 25 ≥ 20, the +first relaxed attempt) — recovers the real empirical PSF (correct FWHM) +instead of falling back, while `relaxation_attempts=0` on the same data +still fails cleanly, proving the relaxation is what does it. Writing that +test surfaced a second thing worth knowing: `fit_moffat_psf`'s own stamp +extraction (`stamp_size=25`, so a ±12 px half-window) doesn't check for +neighbor contamination the way `estimate_psf`'s isolation filter does — +stars closer than ~24 px apart corrupt each other's Moffat fit (tested +directly: clustering the existing fallback test's synthetic stars into a +6 px box made every fit fail to converge, where the original ~25 px +spacing fits cleanly) — not fixed here, since `fit_moffat_psf`'s own +docstring already documents that it deliberately skips the isolation +check as a defensible trade for a *fit* (a bad stamp shows up as a poor +residual, in principle), but the synthetic evidence says that trade has +a real limit worth knowing about. + +**`zogy_subtract`'s `V_ast` opt-in was silent — a direct caller who +simply didn't pass `n_sources`/`r_sources` got `V_ast = 0` with no +signal anything was skipped**, unlike `run_pipeline`, which always +supplies them. The design itself (opt-in at the low-level function, +mandatory at the high-level one) is still the right call — computing +`n_sources`/`r_sources` costs two extra `detect_sources` passes, real +cost a direct caller might legitimately not want. What was missing was +visibility: added a `@warn` (once per session, via `maxlog=1`, so it +doesn't spam a caller who's already made an informed choice) when both +are left `nothing`. Regression test confirms it fires when omitted and +stays silent when both are supplied. blob - f3561f3dbfcf6336eacebefb6260a842d0908dbc blob + 1ecd945e5d1dba046479a0433e2bcd3bb874dbd3 --- src/psf.jl +++ src/psf.jl @@ -22,56 +22,76 @@ An empirical PSF is used rather than an analytic model Moffat) because ZOGY needs the real PSF shape, wings included, and this keeps `estimate_psf` usable on any survey's data without fitting a model family per instrument. But a field sparser or more crowded than expected -can leave zero stars passing the isolation/saturation filter above; when -that happens and `fallback` is true (the default), a warning is emitted -and [`fit_moffat_psf`](@ref) is used instead — an analytic fit trades the -real PSF's exact shape for something usable at all. `fallback=false` -keeps the hard failure instead. +can leave zero stars passing the isolation/saturation filter at the given +`min_separation`; when that happens, `min_separation` is halved and the +search retried, up to `relaxation_attempts` times (default 2, i.e. +`min_separation`, then `/2`, then `/4`) before giving up on the empirical +approach — a field with a few usable stars at a tighter isolation radius +still gives the *real* PSF shape, which a fallback to an analytic model +never can, so this is tried first rather than jumping straight to it. If +even the most relaxed attempt finds nothing and `fallback` is true (the +default), a warning is emitted and [`fit_moffat_psf`](@ref) is used +instead — an analytic fit trades the real PSF's exact shape for something +usable at all. `fallback=false` keeps the hard failure instead. """ function estimate_psf(image::AbstractMatrix{<:Real}; stamp_size::Integer=25, threshold::Real=20.0, min_separation::Real=40.0, - saturate::Real=Inf, fallback::Bool=true) + saturate::Real=Inf, fallback::Bool=true, + relaxation_attempts::Integer=2) isodd(stamp_size) || throw(ArgumentError("stamp_size must be odd")) half = stamp_size ÷ 2 sources = detect_sources(permutedims(image); threshold=threshold) nx, ny = size(image) + function collect_stamps(sep::Real) + stamps = Matrix{Float64}[] + for s in sources + x, y = round(Int, s.x), round(Int, s.y) + (half < x <= nx - half && half < y <= ny - half) || continue + + isolated = all(hypot(s.x - o.x, s.y - o.y) >= sep + for o in sources if !(o.x == s.x && o.y == s.y)) + isolated || continue + + stamp = image[x-half:x+half, y-half:y+half] + maximum(stamp) < saturate || continue + + background = median(vcat(stamp[1, :], stamp[end, :], stamp[:, 1], stamp[:, end])) + stamp = stamp .- background + + total = sum(stamp) + total > 0 || continue + xs = Float64.(-half:half) + cx = sum(xs .* sum(stamp, dims=2)[:]) / total + cy = sum(xs .* sum(stamp, dims=1)[:]) / total + + centered = _shift_bilinear(stamp, -cx, -cy) + s2 = sum(centered) + s2 > 0 || continue + push!(stamps, centered ./ s2) + end + return stamps + end + stamps = Matrix{Float64}[] - for s in sources - x, y = round(Int, s.x), round(Int, s.y) - (half < x <= nx - half && half < y <= ny - half) || continue - - isolated = all(hypot(s.x - o.x, s.y - o.y) >= min_separation - for o in sources if !(o.x == s.x && o.y == s.y)) - isolated || continue - - stamp = image[x-half:x+half, y-half:y+half] - maximum(stamp) < saturate || continue - - background = median(vcat(stamp[1, :], stamp[end, :], stamp[:, 1], stamp[:, end])) - stamp = stamp .- background - - total = sum(stamp) - total > 0 || continue - xs = Float64.(-half:half) - cx = sum(xs .* sum(stamp, dims=2)[:]) / total - cy = sum(xs .* sum(stamp, dims=1)[:]) / total - - centered = _shift_bilinear(stamp, -cx, -cy) - s2 = sum(centered) - s2 > 0 || continue - push!(stamps, centered ./ s2) + sep = min_separation + for attempt in 0:relaxation_attempts + stamps = collect_stamps(sep) + isempty(stamps) || break + attempt == relaxation_attempts && break + sep /= 2 end if isempty(stamps) if fallback - @warn "no isolated, unsaturated stars found for empirical PSF estimation; " * - "falling back to an analytic Moffat fit" threshold min_separation + @warn "no isolated, unsaturated stars found for empirical PSF estimation, " * + "even after relaxing min_separation; falling back to an analytic " * + "Moffat fit" threshold min_separation relaxation_attempts return fit_moffat_psf(image; stamp_size, threshold, saturate) end - error("no isolated, unsaturated stars found for PSF estimation " * - "(try lowering threshold or min_separation)") + error("no isolated, unsaturated stars found for PSF estimation, even after " * + "relaxing min_separation (try lowering threshold or raising relaxation_attempts)") end combined = dropdims(median(cat(stamps...; dims=3); dims=3); dims=3) blob - 6b49312add8f031c6f09447887419699bab828ef blob + 8564ff052f8902fdd43e741c6930005eb7548059 --- src/zogy.jl +++ src/zogy.jl @@ -53,6 +53,13 @@ and `r_sources` (tables with `x`, `y` columns, e.g. fr if reprojection onto a common grid has already removed essentially all registration error — reasonable here since `build_reference` already reprojects reference frames onto the exact science-frame grid. +[`run_pipeline`](@ref) always supplies both, so this only affects a +direct `zogy_subtract` call: omitting them silently (rather than as a +deliberate, informed choice) risks reading `S_corr`'s significance at +face value near bright stars where real sub-pixel misregistration is +exactly what `V_ast` accounts for — a `@warn` (once per session) fires +when both are left unset, so this stays visible rather than a silent +default. """ function zogy_subtract(n_image::AbstractMatrix{<:Real}, r_image::AbstractMatrix{<:Real}; psf_n::AbstractMatrix{<:Real}, psf_r::AbstractMatrix{<:Real}, @@ -93,6 +100,9 @@ function zogy_subtract(n_image::AbstractMatrix{<:Real} V_r = _variance_map(r_sub, sigma_r, gain_r) var_S = _convolve_variance(V_n, k_n) .+ _convolve_variance(V_r, k_r) + if n_sources === nothing && r_sources === nothing + @warn "n_sources/r_sources not given; V_ast (astrometric registration noise) omitted from S_corr — significance may be overstated near bright stars with real sub-pixel misregistration" maxlog=1 + end if n_sources !== nothing && r_sources !== nothing sigma_x, sigma_y = _astrometric_scatter(n_sources, r_sources) if sigma_x > 0 || sigma_y > 0 blob - eceb44a816f34f836fdba4e35ed55e737ccbd384 blob + 0cd7b1fe4c1ac6031cea977865aefbb310a0d534 --- test/runtests.jl +++ test/runtests.jl @@ -300,16 +300,25 @@ using Reproject true_alpha, true_beta = 3.0, 2.5 # centers all within ~25 px of each other, closer than the default # min_separation=40 — every star fails estimate_psf's isolation - # filter, forcing the analytic fallback. + # filter, forcing the analytic fallback. relaxation_attempts=0 + # below keeps this true regardless of the relaxation added for + # sparser fields (see the "min_separation relaxation" testset) — + # this test is specifically about the fallback itself, not about + # how hard estimate_psf tries before reaching it, and these + # centers are also close enough that a relaxed min_separation + # would make some pass isolation while still leaving overlapping + # fit_moffat_psf stamps (stamp_size=25 => half=12, closer than + # 2*12=24 apart) that would corrupt the analytic fit anyway. centers = [(50, 50), (70, 55), (55, 70), (75, 75)] for (cx, cy) in centers, i in 1:nx, j in 1:ny r2 = (i - cx)^2 + (j - cy)^2 img[i, j] += 3000.0 / (1 + r2 / true_alpha^2)^true_beta end - @test_throws ErrorException estimate_psf(img; threshold=15.0, fallback=false) + @test_throws ErrorException estimate_psf(img; threshold=15.0, fallback=false, + relaxation_attempts=0) - psf = estimate_psf(img; threshold=15.0) # fallback=true by default + psf = estimate_psf(img; threshold=15.0, relaxation_attempts=0) # fallback=true by default @test sum(psf) ≈ 1.0 atol=1e-6 c = size(psf, 1) ÷ 2 + 1 @@ -325,6 +334,40 @@ using Reproject @test_throws ErrorException fit_moffat_psf(empty_image; threshold=15.0) end + @testset "estimate_psf min_separation relaxation" begin + Random.seed!(14) + nx, ny = 120, 120 + img = 100.0 .+ 3.0 .* randn(nx, ny) + true_sigma = 1.8 + # 25 px apart: fails the default min_separation=40, but 25 >= 20 + # (the first relaxed attempt, 40/2), so both should pass isolation + # once relaxed — this should recover a real empirical PSF, not + # fall back to the analytic Moffat fit. + centers = [(40, 60), (65, 60)] + for (cx, cy) in centers, i in 1:nx, j in 1:ny + r2 = (i - cx)^2 + (j - cy)^2 + img[i, j] += 4000.0 * exp(-r2 / (2 * true_sigma^2)) + end + + # relaxation_attempts=0 must not relax at all: with fallback=false + # (so this can't succeed via fit_moffat_psf either), it has to + # fail cleanly at the original min_separation=40. + @test_throws ErrorException estimate_psf(img; threshold=15.0, fallback=false, + min_separation=40.0, relaxation_attempts=0) + + # With relaxation (the default, relaxation_attempts=2), the same + # data succeeds via the empirical path instead of falling back — + # verified by the recovered FWHM matching the injected Gaussian's. + psf_relaxed = estimate_psf(img; threshold=15.0, min_separation=40.0) + @test sum(psf_relaxed) ≈ 1.0 atol=1e-9 + c = size(psf_relaxed, 1) ÷ 2 + 1 + @test argmax(psf_relaxed) == CartesianIndex(c, c) + half_max_pos = findfirst(r -> psf_relaxed[c+r, c] < psf_relaxed[c, c] / 2, 0:(c - 1)) + @test half_max_pos !== nothing + half_max_r = half_max_pos - 1 + @test half_max_r - 1 <= true_sigma * 2sqrt(2log(2)) / 2 <= half_max_r + 1 + end + @testset "zogy_subtract" begin psf = zeros(9, 9) for i in 1:9, j in 1:9 @@ -401,6 +444,16 @@ using Reproject sigma_n=sigma, sigma_r=sigma, gain_n=1e8, gain_r=1e8) @test isapprox(mean(s_corr_bg), 0.0; atol=0.2) @test isapprox(std(s_corr_bg), 1.0; atol=0.15) + + # regression test: omitting n_sources/r_sources (silently dropping + # V_ast) used to be undiscoverable without reading the docstring — + # now warns explicitly, one call with, one without. + sources_n = Table(x=[64.0], y=[64.0]) + sources_r = Table(x=[64.2], y=[64.1]) + @test_logs (:warn, r"V_ast") match_mode=:any zogy_subtract( + img, img; psf_n=psf, psf_r=psf, sigma_n=1.0, sigma_r=1.0) + @test_logs zogy_subtract(img, img; psf_n=psf, psf_r=psf, sigma_n=1.0, sigma_r=1.0, + n_sources=sources_n, r_sources=sources_r) end @testset "run_pipeline with reference" begin