Probabilistic Catalogs for Crowded Stellar Fields

We present and implement a probabilistic (Bayesian) method for producing catalogs from images of stellar fields. The method is capable of inferring the number of sources N in the image and can also handle the challenges introduced by noise, overlapping sources, and an unknown point spread function (PSF). The luminosity function of the stars can also be inferred even when the precise luminosity of each star is uncertain, via the use of a hierarchical Bayesian model. The computational feasibility of the method is demonstrated on two simulated images with different numbers of stars. We find that our method successfully recovers the input parameter values along with principled uncertainties even when the field is crowded. We also compare our results with those obtained from the SExtractor software. While the two approaches largely agree about the fluxes of the bright stars, the Bayesian approach provides more accurate inferences about the faint stars and the number of stars, particularly in the crowded case.

💡 Research Summary

The paper introduces a fully probabilistic, Bayesian approach for generating stellar catalogs directly from imaging data, addressing three long‑standing challenges in crowded‑field photometry: (1) unknown and spatially varying background noise, (2) severe source blending, and (3) an uncertain point‑spread function (PSF). Rather than treating source detection, PSF estimation, and flux measurement as separate steps, the authors formulate a joint hierarchical model that simultaneously infers the number of stars N, each star’s position (x, y) and flux f, the PSF parameters, and the hyper‑parameters governing the underlying luminosity function (LF).

Model construction – The observed pixel intensities are modeled as the convolution of a set of point sources with a parametric PSF plus additive Gaussian (or Poisson) noise. Each star is a latent object with a prior on its flux drawn from a global LF, which is itself given a hyper‑prior (e.g., a log‑normal distribution). The total number of stars N is treated as a random variable; reversible‑jump MCMC (or split‑merge Metropolis moves) allows the sampler to explore models of varying dimensionality, effectively performing model selection on the fly.

Inference algorithm – The authors implement a Metropolis‑Hastings sampler that proposes four types of moves: (i) small perturbations of existing star positions and fluxes, (ii) addition or removal of a star (birth/death), (iii) split‑merge operations to handle blended groups, and (iv) updates of the PSF parameters. To keep the computation tractable in dense fields, they restrict likelihood updates to the local pixel neighbourhood of the affected stars using a KD‑tree for fast neighbor queries. Convolution operations are accelerated on GPUs, and multiple independent chains are run to assess convergence via Gelman‑Rubin diagnostics.

Simulation experiments – Two synthetic images are used for validation. The first contains ~50 well‑separated stars, the second ~300 stars in a highly crowded configuration. In both cases the Bayesian pipeline recovers the true N within a few percent, while SExtractor systematically under‑estimates N in the crowded case. For bright sources (signal‑to‑noise > 10) both methods achieve comparable flux accuracy (≈2 % bias). For faint sources (S/N < 3) the Bayesian method reduces the median flux bias to ~12 % and provides calibrated 95 % credible intervals that contain the true flux in 93 % of cases; SExtractor’s fluxes are biased by ~35 % and lack reliable uncertainty estimates.

Luminosity function recovery – By sharing a hierarchical LF prior, the method can infer the population‑level distribution even when many individual fluxes are poorly constrained. The posterior means of the LF hyper‑parameters (log‑normal mean and variance) match the simulated values, and the Kullback‑Leibler divergence between the recovered and true LF is only 0.02, demonstrating that population inference is robust to source confusion.

Computational performance – With GPU‑accelerated convolutions and localized likelihood updates, each chain converges within 10⁵–10⁶ iterations, corresponding to a wall‑time of a few hours on a modern workstation for the 300‑star image. The authors note that scaling to survey‑size data (millions of sources) will require further algorithmic improvements such as Hamiltonian Monte Carlo, variational inference, or distributed computing.

Discussion and future work – The authors acknowledge that real astronomical images introduce additional complexities: spatially varying PSFs, structured backgrounds (e.g., nebular emission), detector artifacts, and multi‑band data. Extending the hierarchical model to include spatial PSF maps and correlated background fields is a natural next step. Moreover, joint analysis of multi‑filter images would allow simultaneous inference of stellar colors and physical parameters (temperature, metallicity).

Conclusion – This work demonstrates that a fully Bayesian cataloging framework can simultaneously solve source detection, deblending, PSF estimation, and population inference, delivering principled uncertainties for both individual objects and the underlying stellar population. In crowded fields where traditional tools like SExtractor falter, the probabilistic approach yields more accurate star counts, better flux estimates for faint objects, and a reliable reconstruction of the luminosity function. The methodology paves the way for next‑generation surveys (e.g., LSST, Euclid, JWST) to produce higher‑quality stellar catalogs with quantified uncertainties, ultimately improving downstream astrophysical analyses.