commit - 8ff5933c85c68fcb18d5be2ee05c05d4d2c6df00
commit + 92ed176ba250755bdfec88536921f1d0f69f0113
blob - 413f5933626675d76786f401e56e3c9361681dbf
blob + dd66e303cdccc4f3354ebc934d1667928a28cfc4
--- docs/src/index.md
+++ docs/src/index.md
## 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
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
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
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},
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
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
@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
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
