PaperlayerView exact revision on arXiv

Source-linked AI summary

GalSim: The modular galaxy image simulation toolkit

Barnaby Rowe, Mike Jarvis, Rachel Mandelbaum, Gary M. Bernstein, James Bosch, Melanie Simet, Joshua E. Meyers, Tomasz Kacprzak, Reiko Nakajima, Joe Zuntz, Hironao Miyatake, Joerg P. Dietrich, Robert Armstrong, Peter Melchior, Mandeep S. S. Gill

arXiv:1407.7676v3astro-ph.IMastro-ph.CO

TL;DR

Astronomical image analysis needs simulations because weak-lensing and related measurements can contain systematic biases requiring high-precision control. GALSIM provides a modular open-source toolkit for generating, transforming, rendering, and analyzing simulated images, and tests show Stage IV-level accuracy across several demanding cases. It can serve as a benchmark and extensible reference codebase, but users should test it for their specific applications.

  • Problem

    Astronomical image analysis needs simulations because weak-lensing and related measurements can contain systematic biases requiring high-precision control.

  • Method

    GALSIM provides a modular open-source toolkit combining astronomical object models, transformations, rendering, noise models, and analysis tools.

  • Results

    GALSIM output meets Stage IV weak-lensing accuracy requirements for analytic Sérsic, image-derived, and reconvolution cases.

  • Takeaways & Limitations

    GALSIM can serve as an open-source benchmark, common reference point, and extensible codebase for astronomical image-analysis development.

  • Takeaways & Limitations

    The reported tests cannot validate GALSIM for every future application, so users should conduct tests tailored to their requirements.

Abstract

from arXiv · show

GALSIM is a collaborative, open-source project aimed at providing an image simulation tool of enduring benefit to the astronomical community. It provides a software library for generating images of astronomical objects such as stars and galaxies in a variety of ways, efficiently handling image transformations and operations such as convolution and rendering at high precision. We describe the GALSIM software and its capabilities, including necessary theoretical background. We demonstrate that the performance of GALSIM meets the stringent requirements of high precision image analysis applications such as weak gravitational lensing, for current datasets and for the Stage IV dark energy surveys of the Large Synoptic Survey Telescope, ESA's Euclid mission, and NASA's WFIRST-AFTA mission. The GALSIM project repository is public and includes the full code history, all open and closed issues, installation instructions, documentation, and wiki pages (including a Frequently Asked Questions section). The GALSIM repository can be found at https://github.com/GalSim-developers/GalSim .

1. Introduction

Increasing survey data volumes reduce statistical uncertainty while demanding tighter control of systematic errors in astronomical image analysis. GALSIM was conceived to simulate these effects, especially for weak-lensing measurements and precision cosmology.

  • Larger surveys increase galaxy counts and reduce statistical uncertainty, requiring proportionally smaller systematic errors in inferred galaxy properties.
  • Weak-lensing shear estimates from noisy galaxy images can suffer generic systematic biases unless measurement methods are constructed carefully.
  • Model bias arises when true galaxy surface-brightness profiles differ from the fitted models, motivating tests across shear-measurement methods.
  • Precision photometry and shape estimation also face object confusion, chromatic mismatches, atmospheric refraction, colour gradients, and detector nonlinearities.
  • Simulations are needed to estimate or calibrate shear biases and demonstrate that weak-lensing and related measurements meet precision requirements.

2. Software overview

GALSIM is an open-source, modular class-library toolkit for combining astronomical object models, transformations, lensing, rendering, noise, and shape estimation. Its development emphasizes extensibility, documentation, testing, review, and numerical validation.

  • Capabilities: The toolkit represents astronomical objects, applies transformations and convolutions, generates lensing distortions, maps image and world coordinates, renders images, adds noise, and estimates shapes.
  • Scientific scope: The toolkit targets accurate photometry and astrometry as well as weak-lensing simulations, making it useful for applications such as photometric-redshift estimation.
  • Software design: GALSIM was designed as a modular, extensible toolkit whose components can be maintained independently and used interchangeably through common interfaces.
  • Development process: New features require documentation and example scripts, while new code requires tests and extensive review of code, documentation, and unit tests.
  • Interface: GALSIM is fundamentally a Python class library with numerical calculations implemented largely in C++, while shielding users from the C++ layer.
  • Interface: Users can perform most major operations through standalone executables that read configuration files as an alternative to writing Python code.

3. Surface brightness profiles

GALSIM provides analytic and image-based surface-brightness models for galaxies and stars, alongside optical and atmospheric PSF models. Profiles expose physically motivated parameters and support transformations, truncation, interpolation, and high-precision rendering.

  • Galaxy profiles: GALSIM includes analytic galaxy profiles such as exponential, de Vaucouleurs, and Sérsic models, with sizes expressible using scale or half-light radii.
  • Galaxy profiles: De Vaucouleurs profiles have cuspy cores and broad wings, so GALSIM supports specifying a truncation radius to limit their extent.
  • Galaxy profiles: The Sérsic profile generalizes exponential and de Vaucouleurs profiles, with GALSIM restricting its index to 0.3 ≤ n ≤ 6.2 for numerical and rendering accuracy.
  • Image-based profiles: Observed HST galaxy images can serve as realistic light-profile models through GALSIM’s RealGalaxy class, preserving morphological variation and irregular galaxies.
  • PSF models: Optical PSFs include Airy diffraction, central obscurations, support struts, and Zernike-mode aberrations through the Airy and OpticalPSF classes.
  • PSF models: OpticalPSF profiles are generated by discrete Fourier transforming the PSF into real space and interpolating the resulting image.
  • PSF models: The Moffat model represents stellar PSFs with broader wings than a Gaussian and approaches a Gaussian as β →∞.

4. Lensing shear and magnification

GALSIM applies weak-lensing shear and convergence transformations and can generate coherent fields from power spectra or halo models. Its Fourier-grid implementation includes safeguards for finite sampling, aliasing, and interpolation effects.

  • Lensing transformations: GALSIM models shear as anisotropic stretching and convergence as isotropic changes in apparent object size, with magnification produced by their combination.
  • Lensing transformations: GALSIM also offers an area-conserving shear transformation in addition to the non-area-conserving weak-lensing transformation.
  • Lensing fields: The toolkit supports single shear and convergence values, cosmological power-spectrum fields, spherical NFW-halo fields, and user-supplied E- and B-mode spectra.
  • Lensing fields: The default cosmological field generator uses a Gaussian random-field approximation and is intended for semi-realistic fields rather than highly accurate end-to-end cosmological inference.
  • Fourier implementation: The shear-generation algorithm samples Fourier amplitudes and random phases on a finite grid, then transforms them back to real space with periodic boundary conditions.
  • Numerical limitations: Finite k ranges can omit power or cause aliasing, so GALSIM supports grid expansion, correlation-function prediction, and Nyquist-scale band-limiting.
  • Interpolation limitations: Interpolation modifies shear two-point functions by more than 10% below three original grid spacings, making grid spacing and edge effects important.

5. World coordinate systems

GALSIM maps astronomical profiles between sky and image coordinates through flexible WCS models, using local affine transformations and flux-preserving pixel integration during rendering.

  • Coordinate-system types: GALSIM supports WCS relationships ranging from a simple pixel scale to complex functions like those recorded in real FITS headers.Celestial systems use RA and Dec; Euclidean systems use tangent-plane coordinates (u, v).
  • Coordinate-system types: The toolkit represents celestial WCS with user-defined or FITS-derived mappings and Euclidean WCS with arbitrary functions or affine transformations.PixelScale is the simplest affine specialization and describes uniform square pixels.
  • Local transformations: GALSIM converts profiles locally by tangent projection and Jacobian evaluation, assuming WCS derivatives are approximately constant across each object.The Jacobian defines the local affine approximation used over the profile’s extent.
  • Local transformations: Transforming from world to image coordinates applies J^-1 and multiplies the profile by |det J| to preserve total flux.The toolkit’s toImage and toWorld functions perform the required conversions in either direction.
  • Pixel and PSF treatment: Realistic rendering requires consistent coordinate choices for galaxy profiles, pixel responses, and PSFs, including avoiding a second pixel convolution for measured PSFs.The no_pixel option supports rendering without the extra pixel convolution.
  • Pixel integration: Pixel integration convolves the profile with the pixel response at each pixel centre and multiplies by pixel area, a calculation valid in either coordinate system.For complicated WCS functions, drawImage handles the coordinate conversion and convolution automatically.

6. Image rendering

GALSIM offers direct, Fourier, and photon-shooting rendering methods, selecting among them according to object properties and desired pixel treatment while documenting method-specific approximations and limits.

  • Rendering methods: GALSIM renders objects directly in real space, through discrete Fourier transforms, or by photon shooting from the surface-brightness distribution.The methods are theoretically equivalent apart from photon-shot noise, with additional differences from technique-specific approximations.
  • Rendering methods: Photon shooting automatically integrates profiles over pixels by binning sampled photons, whereas no_pixel uses direct rendering or DFT.The photon method treats the surface brightness as a probability distribution and samples a finite number of photons.
  • Method limitations: Rendering capabilities differ by method: real-space convolution is restricted, deconvolution requires DFT, and flexion cannot currently be rendered in Fourier space.Direct rendering can handle convolution of only two profiles, while some photon-shooting cases are not implemented.
  • Pixel treatment: The default drawImage operation convolves an object with a square pixel response before normally rendering the resulting convolution with DFT.Direct centre sampling is available but does not correspond directly to a real image and may not conserve input flux in the pixel sum.
  • Interpolated images: InterpolatedImage direct rendering is usually slow because each output pixel sums over multiple input pixels according to the chosen interpolant.Interpolation extends the object footprint and changes the effective Fourier cutoff; Lanczos interpolants offer more compact support than sinc.
  • Photon shooting: Photon-shooting accuracy depends on narrow annuli, and GALSIM reports close correspondence between photon-shooting and DFT results in numerical validation.The procedure uses Nγ photon locations and weights, with additional evaluations determined by the rejected-photon fraction ε.
  • Realistic galaxies: Reconvolution removes the HST PSF, applies shear or magnification, and renders high-resolution real galaxies as they would appear through a lower-resolution telescope.GALSIM also models correlated noise through reconvolution and whitens the final image to uncorrelated Gaussian noise.

7. Noise models

GALSIM provides stochastic models for uncorrelated, nonstationary, detector-like, and correlated image noise. Its COSMOS noise model and whitening procedure address inter-pixel correlations while documenting variance and signal-to-noise trade-offs.

  • GALSIM provides stochastic noise models because astronomical images contain photon, thermal, electronic, and other noise sources.
  • GaussianNoise adds stationary, uncorrelated N(0, σ2) noise, while PoissonNoise varies its variance with each pixel’s expected counts λi.
  • CCDNoise combines Poisson photon-count noise with stationary Gaussian read noise and optionally converts photon counts to ADU using a gain.
  • GALSIM models correlated noise through a discrete autocorrelation function whose stationary form determines the covariance matrix for all pixel pairs.The function is indexed by integer pixel separations n and m and is peaked at zero separation.
  • Noise whitening adds a designed noise field so the combined image noise becomes approximately uncorrelated and stationary, with COSMOS support through getCOS-MOSNoise().The supplied COSMOS model represents correlations measured from blank sky regions in COSMOS F814W images.
  • ∼10−2–10−3 of the zero-lag variance: whitening reduced interpixel covariances to this level in a ∼400 × 300 pixel COSMOS patch.For sheared and convolved COSMOS fields, correlations were reduced by 2–3 orders of magnitude and became statistically undetectable except in the largest simulations.
  • Whitening increases output variance, while its information cost depends on the measured property, object profile, and original noise correlations.For transformed fields, the added whitening variance depends on shear, PSF, and final pixel scale and is typically lower than for raw native-resolution COSMOS noise.

8. Shear estimation

GALSIM includes image-based routines for weighted galaxy moments and PSF-corrected shear estimation. These routines implement established methods and add interface, computational, masking, and correctness improvements.

  • GALSIM includes routines for weighted galaxy-shape moments and PSF-corrected shear estimation, including the re-Gaussianization method.
  • FindAdaptiveMoments iteratively fits an elliptical Gaussian and uses that fitted profile as the weight function for adaptive moments.
  • EstimateShear implements multiple PSF-correction methods, including kurtosis-based, linear non-Gaussianity, and re-Gaussianization approaches.
  • The first three PSF-correction methods output per-galaxy distortion, which requires an ensemble responsivity factor to obtain shear estimates.Their outputs are explicitly distinct from per-galaxy shear estimates produced by KSB.
  • Users can tune overall and method-specific behavior through HSMParams and arguments passed to individual functions.
  • The GALSIM implementations improve usability and speed, optimize convolutions and transforms, support masks and extensible weight maps, and fix nonsquare-input-PSF handling.

9. Numerical validation

GALSIM validation uses adaptive-moment measurements of rendered ellipticity and size to test whether image-rendering discrepancies remain below stringent weak-lensing requirements. Across DFT rendering, interpolated-image transformations, and reconvolution, the tested defaults generally satisfy these targets, while limitations require application-specific validation.

  • Validation framework: The validation metric measures rendering accuracy through adaptive-moment estimates of object ellipticity and size.Errors are modeled for ellipticity components g1 and g2 and size σ.
  • Validation framework: 2 × 10^-4 and 2 × 10^-5 define the target limits for multiplicative and additive discrepancies, respectively, at one tenth of Euclid requirements.The multiplicative target applies to m-type biases, while the additive target applies to c-type biases.
  • DFT rendering: DFT rendering met the m-type bias tolerances across the tested Sérsic indices, while c-type additive biases were consistent with zero.The comparison used photon-shooting and DFT-rendered Sérsic profiles with differences evaluated in ellipticity and size.
  • Interpolated images: Setting pad_factor = 4 and using a two-dimensional quintic interpolant satisfied mi,interp < 2 × 10^-4 and ci,interp < 2 × 10^-5 across the tested transformations.Shear-component leakage was also consistent with zero and comfortably within requirements for pad_factor values of 4 and 6.
  • Reconvolution: At the default folding_threshold of 5 × 10^-3, reconvolution met the multiplicative requirement in ensemble averages, but 3/200 galaxies exceeded the m2,reconv limit.Reducing folding_threshold by a factor of 10 brought m2,reconv within requirements for all 200 tested galaxies; additive biases remained extremely small.
  • Limitations: The validation suite cannot establish performance for all future applications, so users should test configurations tailored to their data and requirements.The tests used circular PSFs and a particular COSMOS training sample, limiting coverage of nonsymmetric PSFs and other galaxy populations or filters.

10. Performance

GALSIM prioritizes accurate rendering while using C++ acceleration, specialized transforms, adaptive integration, FFT optimizations, and configurable parallelism to achieve practical simulation speeds. Typical analytic profiles render in about 0.01 seconds per object, while interpolated images take about 0.1 seconds and photon shooting can be faster for faint objects.

  • Optimizations: GALSIM prioritizes rendering accuracy and accelerates time-critical operations through C++ called from its Python interface.The implementation uses standard compiled-code optimizations, including precomputation, tight loops, vectorization-friendly iterators, lookup tables, and caches.
  • Optimizations: Radially symmetric profiles can use one-dimensional Hankel transforms instead of two-dimensional FFTs, with particular speed benefits for high-n Sérsic profiles.The transform is precomputed in one dimension and used to fill Fourier-space values.
  • Optimizations: Adaptive Gauss-Kronrod-Patterson integration supports Hankel transforms, real-space integration, photon-shooting probabilities, and Sérsic half-light-radius calculations.The algorithm adds points until convergence or 175 points, then recursively subdivides unconverged regions.
  • Optimizations: Accurate Sine integral evaluation was developed for Lanczos-interpolated images because standard formulae reached only about 10^-6 accuracy.GALSIM’s formulae achieve 10^-16 accuracy and are given in Appendix B.
  • Timings: Most analytic profiles render in about 0.01 seconds per object, whereas interpolated-image models typically take about 0.1 seconds per object.Interpolated-image timing depends strongly on the defining image size.
  • Timings: Photon shooting scales linearly with photon count and is typically faster than DFT rendering below a crossover of 10^3–10^4 photons.The crossover depends on the profile, and photon shooting can be significantly faster for faint objects.
  • Pitfalls: RealGalaxy whitening can increase rendering time to about 1 second per object, depending on pixel scale and PSF.Whitening removes correlated noise from the original HST images after PSF convolution.

11. Effects not in GALSIM

GALSIM does not yet model several observational and instrumental effects that may matter for future extragalactic surveys. These gaps include atmospheric, detector, artifact, contaminant, and non-affine lensing effects.

  • Unmodeled effects: GALSIM lacks physical multilayer atmospheric PSF models, despite their potential use in future survey simulations.Adding such a module was investigated during GREAT3 preparation, but the functionality remains absent.
  • Unmodeled effects: GALSIM cannot model detector saturation, bleeding, interpixel capacitance, or other cross-talk effects that affect image flux and shape determination.These are grouped among nonlinear detector effects not yet supported by the toolkit.
  • Unmodeled effects: GALSIM does not simulate ghosts, detector persistence, cosmic rays, satellite trails, or meteors.The paper identifies these artifacts and near-field contaminants as relevant to reliable object measurements and upcoming surveys.
  • Unmodeled effects: Flexion is not implemented because non-affine transformations are incompatible with GALSIM’s Fourier-space rendering of individual profiles.The standalone GALFLEX module can be integrated through a NumPy array interface.
  • Unmodeled effects: GALSIM currently omits spatial variation in photon-detection efficiency, including intra-pixel quantum-efficiency variation.Large-scale effects could be addressed through post-processing, while significant intra-pixel variation may require careful rendering design.
  • Unmodeled effects: Complicated detector astrometric shifts such as edge distortions and tree-rings are not modeled by existing WCS packages and would require custom implementation.These effects are listed as examples of complicated WCS functions absent from GALSIM.

12. Conclusions

The paper presents GALSIM as an open, modular toolkit for high-precision astronomical image simulation, combining flexible physical models, multiple rendering strategies, and noise handling. Tests show Stage IV-level accuracy for weak-lensing transformations, while users must tailor validation and account for effects GALSIM does not yet simulate.

  • Conclusions: GALSIM is an open, transparent, community-maintained toolkit for high-precision simulations of deep extragalactic imaging surveys.The paper describes its physical models, world-coordinate support, rendering strategies, and noise addition.
  • Conclusions: GALSIM brings novel photon shooting, noise whitening, and consistent transformations under varied WCS transformations together in a flexible class-library interface.The paper identifies this combination as unique within astronomical image simulation.
  • Conclusions: Multiple parallel rendering methods were important for testing and validating GALSIM output.The paper describes this capability as unique among comparable simulation packages.
  • Conclusions: GALSIM tests demonstrate weak-lensing shear and magnification accuracy at the level required for Euclid, LSST, and WFIRST-AFTA.The demonstrations cover Sérsic, InterpolatedImage, RealGalaxy, and reconvolution cases.
  • Conclusions: Rendering accuracy can be traded against execution speed through GSParams, so users should modify validation tests for application-specific requirements.The paper notes that no validation suite can cover every usage scenario.
  • Conclusions: Important effects remain unsimulated and may significantly affect data analysis for upcoming surveys, with continued development expected.The paper frames GALSIM as an ongoing project supported by growing numbers of users and developers.
  • Conclusions: GALSIM can serve as a benchmark, code-comparison reference, and extensible codebase for astronomical image-analysis applications.The paper reports that it has already supported GREAT3 image generation and meets demanding rendering tolerances.

Appendix A.1. Grid parameters

Appendix A.1 defines the finite real-space and Fourier-space grids used to approximate continuous functions with discrete Fourier methods. Grid length and spacing determine the number of samples and the sampled Fourier range under the paper’s angular-frequency convention.

  • Grid parameters: GALSIM approximates continuous functions on a finite, discrete real-space grid characterized by its length, spacing, and number of points.The number of grid points along one dimension is N = L/d.
  • Grid parameters: The Fourier-space grid has the same N as the real-space grid, with spacing kmin = 2π/L.This spacing corresponds to a Fourier mode that exactly fits inside the square grid.
  • Grid parameters: The one-dimensional Fourier range extends from k1 = −π/d to π/d, corresponding to a maximum sampled wavenumber kmax = π/d.The upper limit represents a mode sampled exactly twice.
  • Transform convention: The paper uses a non-unitary Fourier-transform convention with angular rather than normal frequency in the exponent.This convention introduces additional factors of 2π compared with common signal-processing and DFT conventions.
  • Transform convention: The appendix introduces Poisson summation as the key identity for approximating continuous functions with discrete Fourier transforms.The stated form follows the paper’s Fourier-transform conventions.

Appendix A.3. Fourier transform of discrete samples of the power spectrum

The appendix develops a discrete Fourier representation for sampled power-spectrum functions, using periodicity, sampling, truncation, and inverse-transform conventions to support numerical calculation.

  • Approximation: For sufficiently small ∆k, periodic copies of the correlation function are widely separated, allowing non-overlapping regions to approximate ξ+(θ1, θ2).The periodic summation has period L in the real domain.
  • Discrete transforms: The resulting discrete relations are inverse discrete Fourier transforms under the paper’s non-unitary convention.The appendix also relates these expressions to standard inverse DFT and DTFT conventions.
  • Periodic representation: The power-spectrum expression is periodic with period L = 2π/∆k, so one period contains all information about ξ+(θ1, θ2).This periodicity motivates sampling over a single real-space period.
  • Periodic representation: N equally spaced samples over one period have real-space spacing ∆θ = 2π/∆kN = d in each dimension.The sampled function is defined by substituting these grid locations into the periodic Fourier expression.
  • Discrete transforms: Practical computation truncates the infinite Fourier-space sequence to finite length, commonly using N samples in both real and Fourier space.This truncation defines the finite sampled power-spectrum array used numerically.

Appendix A.4. Generating Gaussian fields

The appendix generates Gaussian random lensing fields by sampling convergence modes in Fourier space, transforming them into shear, and applying an inverse DFT to obtain a real-space grid realization.

  • Fourier-space generation: A complex, 2D Fourier-space array is constructed to realize a Gaussian random lensing field with ensemble-average power spectrum approximately P(k).The construction begins from the desired Fourier-space power-spectrum representation.
  • Fourier-space generation: Each Fourier-space convergence mode uses independent standard Gaussian random deviates at every array location.The random field is complex-valued, with real and imaginary Gaussian components.
  • Shear construction: The complex convergence field is converted into a Fourier-space shear field using the phase of the wavevector k = ∆k(p + iq).The E- and B-mode convergence fields are treated analogously.
  • Real-space realization: Applying the inverse DFT to the Fourier-space shear produces a real-space grid g[n, m], whose real and imaginary parts are the two shear components.The resulting field is correctly scaled to have discrete correlation function ξ∆θ[n, m].
  • Real-space realization: The construction yields a field with the target discrete correlation function by design.This follows from the inverse-DFT relation used to map Fourier-space modes to the real-space grid.

Appendix B. Efficient evaluation of the Sine and Cosine integrals

The appendix develops efficient numerical evaluations of the sine and cosine integrals used in GALSIM, combining small- and large-argument approximations to achieve high accuracy and efficiency.

  • Purpose: GALSIM requires efficient evaluation of Si(x) for calculations involving Lanczos interpolants, while Ci(x) is provided as an accurate completeness result.The authors developed these calculations because they could not find a suitably accurate efficient source for Si(x).
  • Small arguments: For small arguments, Padé approximants are applied to the convergent Taylor series of the trigonometric integrals.The resulting formulae are accurate to better than 10−16 for 0 ≤ x ≤ 4.
  • Small arguments: The small-argument formulae achieve accuracy better than 10−16 over 0 ≤ x ≤ 4.Maple was used to derive the displayed approximations.
  • Large arguments: For large arguments, helper functions based on asymptotic behavior are used to construct efficient approximants.The relevant limits are f(x) → 1/x and g(x) → 1/x2 as x →∞.
  • Large arguments: The large-argument Chebyshev-Padé expansions are accurate to better than 10−16 for x ≥4.GALSIM evaluates the resulting polynomials using Horner’s rule for efficiency.
Loading 1407.7676v3…