Image Simulation

STDPipe provides comprehensive image simulation capabilities for generating realistic astronomical images with stars, galaxies, and various imaging artifacts. This is particularly useful for:

• Testing photometry and detection algorithms

• Training real-bogus classifiers for transient detection

• Estimating detection limits and completeness

• Validating pipeline performance

The simulation module provides both high-level and incremental APIs for maximum flexibility.

Overview

The stdpipe.simulation module can generate:

Astronomical Sources:

• Stars with Gaussian or Moffat PSF (pixel-integrated for accuracy)

• Galaxies with Sersic profiles (exponential, de Vaucouleurs, etc.)

Imaging Artifacts:

• Cosmic rays (sharp, tapered, or worm-like tracks)

• Hot pixels (isolated or clustered)

• Bad columns/rows (dead, hot, or noisy)

• Satellite trails (with optional tumbling)

• Diffraction spikes (automatically added to bright stars)

• Optical ghosts (reflections from very bright sources)

All simulated sources are tracked in a catalog with detailed metadata, making it easy to train machine learning classifiers or validate detection algorithms.

Quick Start

The simplest way to create a simulated image is using the high-level simulate_image() function:

from stdpipe import simulation

# Create a realistic simulated image
result = simulation.simulate_image(
    width=500,
    height=500,
    n_stars=100,
    n_galaxies=20,
    n_cosmic_rays=10,
    background=1000.0,
    readnoise=10.0
)

image = result['image']      # The simulated image
catalog = result['catalog']  # Catalog of all injected sources

print(f"Created {len(catalog)} objects")
print(f"Image shape: {image.shape}")

The catalog contains all injected sources with their properties and a type column to distinguish between different object types.

High-Level API

Single-Call Image Simulation

The simulate_image() function creates a complete simulated image in a single call:

from stdpipe import simulation
import numpy as np

# Create a complex simulated field
result = simulation.simulate_image(
    width=1000,
    height=1000,
    # Stars
    n_stars=200,
    star_flux_range=(100, 50000),
    star_fwhm=3.5,
    star_psf='moffat',           # or 'gaussian'
    star_beta=2.5,               # Moffat beta parameter
    # Galaxies
    n_galaxies=50,
    galaxy_flux_range=(500, 10000),
    galaxy_r_eff_range=(3, 15),
    galaxy_n_range=(0.5, 4.0),   # Sersic index
    galaxy_ellipticity_range=(0.0, 0.7),
    # Artifacts
    n_cosmic_rays=15,
    cosmic_ray_length_range=(10, 60),
    n_hot_pixels=30,
    n_satellites=2,
    satellite_length_range=(200, 500),
    n_bad_columns=1,
    # Advanced features
    diffraction_spikes=True,     # Add spikes to bright stars
    spike_threshold=50000,       # Flux threshold for spikes
    optical_ghosts=True,         # Add ghosts to very bright stars
    ghost_threshold=100000,
    # Noise model
    background=1000.0,
    readnoise=10.0,
    gain=1.0,
    # Geometry
    edge=20,                     # Keep sources away from edges
    wcs=None,                    # Optional WCS for RA/Dec
    verbose=True
)

image = result['image']
catalog = result['catalog']
background = result['background']
noise = result['noise']

# Analyze the catalog
stars = catalog[catalog['type'] == 'star']
galaxies = catalog[catalog['type'] == 'galaxy']
artifacts = catalog[catalog['is_real'] == False]

print(f"Stars: {len(stars)}")
print(f"Galaxies: {len(galaxies)}")
print(f"Artifacts: {len(artifacts)}")

The is_real column in the catalog distinguishes astronomical sources (True) from artifacts (False), which is perfect for training real-bogus classifiers.

Incremental API

For more control, you can build images incrementally by adding components one at a time.

Adding Stars

Use add_stars() to add stars with customizable PSF:

from stdpipe import simulation
import numpy as np

# Create blank image
image = np.zeros((500, 500))

# Add stars with Gaussian PSF
catalog_gauss = simulation.add_stars(
    image,
    n=50,
    flux_range=(1000, 20000),
    fwhm=3.0,
    psf_type='gaussian',
    edge=20
)

# Add more stars with Moffat PSF
catalog_moffat = simulation.add_stars(
    image,
    n=30,
    flux_range=(500, 10000),
    fwhm=4.0,
    psf_type='moffat',
    beta=2.5,
    diffraction_spikes=True,    # Add spikes to bright stars
    spike_threshold=15000,
    optical_ghosts=True,        # Add ghosts to very bright stars
    ghost_threshold=50000
)

print(f"Added {len(catalog_gauss)} Gaussian stars")
print(f"Added {len(catalog_moffat)} Moffat stars")

Note

The pixel-integrated PSF models eliminate systematic flux biases that can affect photometry testing. For Gaussian PSF, the bias is reduced from +11.6% at FWHM=1.5 to <0.1% for all FWHM values.

Adding Galaxies

Use add_galaxies() to add extended sources with Sersic profiles:

# Add galaxies with various morphologies
catalog = simulation.add_galaxies(
    image,
    n=20,
    flux_range=(1000, 10000),
    r_eff_range=(5, 20),          # Effective radius in pixels
    n_range=(0.5, 4.0),           # Sersic index
    ellipticity_range=(0.0, 0.8),
    edge=50,                      # Keep away from edges
    gain=1.0
)

# Examine galaxy properties
for galaxy in catalog:
    print(f"Galaxy at ({galaxy['x']:.1f}, {galaxy['y']:.1f}): "
          f"n={galaxy['sersic_n']:.1f}, r_eff={galaxy['r_eff']:.1f}, "
          f"e={galaxy['ellipticity']:.2f}, PA={galaxy['position_angle']:.1f}°")

The Sersic index n controls the profile shape:

• n=0.5: Gaussian-like (very compact)

• n=1.0: Exponential disk (spiral galaxies)

• n=4.0: de Vaucouleurs (elliptical galaxies)

Adding Artifacts

Add various imaging artifacts to test detection and classification algorithms:

# Add cosmic ray tracks
cat_cosmic = simulation.add_cosmic_rays(
    image,
    n_rays=10,
    length_range=(10, 50),
    width_range=(1, 3),
    intensity_range=(5000, 20000),
    profile='sharp'              # or 'tapered', 'worm'
)

# Add hot pixels
cat_hot = simulation.add_hot_pixels(
    image,
    n_pixels=20,
    intensity_range=(1000, 5000),
    clustering=True,             # Create clusters
    cluster_size=3
)

# Add satellite trails
cat_satellites = simulation.add_satellite_trails(
    image,
    n_trails=2,
    length_range=(200, 400),
    width_range=(3, 8),
    intensity_range=(10000, 30000),
    profile='gaussian',          # or 'linear'
    tumbling_prob=0.3            # Fraction with intensity variations
)

# Add bad columns
cat_bad = simulation.add_bad_columns(
    image,
    n_columns=1,
    bad_type='dead',             # or 'hot', 'noisy'
    orientation='vertical'       # or 'horizontal'
)

PSF Models

Gaussian PSF

The default PSF model is an oversampled Gaussian in PSFEx-compatible format:

from stdpipe import simulation, psf

# Create a Gaussian PSF model
psf_model = simulation.create_psf_model(
    fwhm=3.5,
    psf_type='gaussian',
    oversampling=2               # Oversampling for sub-pixel accuracy
)

# Get a pixel-level stamp at a specific position
stamp = psf.get_psf_stamp(psf_model, x=12.3, y=12.7, normalize=True)
print(f"PSF sum: {stamp.sum():.6f}")  # Should be ~1.0

Moffat PSF

The Moffat PSF has broader wings than Gaussian, which is more realistic for many telescopes:

# Create a Moffat PSF model
psf_model = simulation.create_psf_model(
    fwhm=4.0,
    psf_type='moffat',
    beta=2.5                     # Beta parameter (DAOPHOT 'moffat25')
)

The beta parameter controls the wing profile:

• beta=1.5: Very extended wings

• beta=2.5: Standard (DAOPHOT ‘moffat25’)

• beta=4.0: Compact, closer to Gaussian

Optical Aberrations

For realistic testing of photometry algorithms under challenging conditions, simulate PSFs with optical aberrations using Zernike polynomial-based wavefront modeling. This is particularly useful for testing algorithm robustness and training real-bogus classifiers with degraded PSFs.

Basic Usage

from stdpipe import simulation

# Create aberrated PSF
psf_aberrated = simulation.create_psf_model(
    fwhm=3.5,
    psf_type='gaussian',
    defocus=1.0,           # 1 wave of defocus
    astigmatism_x=0.5,     # 0.5 waves at 0°
    wavelength=550e-9,     # 550nm (V-band)
)

# Use with star simulation
image = np.zeros((500, 500))
cat = simulation.add_stars(image, n=50, psf=psf_aberrated)

Aberration Parameters

All aberration coefficients are specified in units of wavelengths (waves):

• defocus (float): Defocus aberration (Zernike Z₄). Default: 0.0

  • Typical moderate: 0.5-2.0 waves

  • Severe: >3.0 waves

  • Creates ring/donut shaped PSF

  • Simulates out-of-focus telescope or atmosphere

• astigmatism_x (float): Astigmatism at 0° orientation (Zernike Z₅). Default: 0.0

• astigmatism_y (float): Astigmatism at 45° orientation (Zernike Z₆). Default: 0.0

  • Typical: 0.3-1.5 waves

  • Creates elliptical PSF

  • Orientation varies with focus position

  • Simulates asymmetric optics or atmospheric turbulence

• coma_x (float): Horizontal coma (Zernike Z₇). Default: 0.0

• coma_y (float): Vertical coma (Zernike Z₈). Default: 0.0

  • Typical: 0.2-1.0 waves

  • Creates comet-like tail asymmetry

  • Common in off-axis telescope configurations

• wavelength (float): Wavelength in meters. Default: 550e-9 (V-band)

  • 450e-9: B-band (blue)

  • 550e-9: V-band (green/visual)

  • 650e-9: R-band (red)

  • Affects diffraction pattern scale

• pupil_diameter (float): Normalized pupil diameter. Default: 1.0

  • Controls relative size of aperture

  • Affects diffraction-limited resolution

Physical Basis

Aberrations are modeled using Zernike polynomials (Noll 1976 ordering) via Fourier optics:

1. Wavefront: W(ρ,θ) = Σ cⱼ·Zⱼ(ρ,θ) where cⱼ are aberration coefficients

2. Complex pupil: P = aperture · exp(2πi·W/λ)

3. PSF intensity: |FFT(P)|²

4. Seeing convolution: PSF ⊗ Gaussian/Moffat (atmospheric effects)

Aberrations combine linearly in the wavefront (not PSF intensity). The final PSF represents both diffraction effects from the aberrated optics and atmospheric seeing (controlled by fwhm).

Example: Testing Photometry Robustness

from stdpipe import simulation, photometry, photometry_measure
import numpy as np

# Test photometry with various aberration levels
defocus_levels = [0.0, 0.5, 1.0, 2.0, 3.0]
flux_errors = []

for defocus in defocus_levels:
    # Create aberrated PSF
    psf_model = simulation.create_psf_model(
        fwhm=4.0,
        psf_type='gaussian',
        defocus=defocus,
        size=51
    )

    # Simulate stars
    result = simulation.simulate_image(
        width=500,
        height=500,
        n_stars=50,
        star_flux_range=(5000, 10000),
        star_fwhm=4.0,
        star_psf=psf_model,
        background=1000.0,
        readnoise=10.0
    )

    image = result['image']
    true_cat = result['catalog']

    # Detect and measure
    detected = photometry.get_objects_sep(image, thresh=5.0)
    measured = photometry_measure.measure_objects(
        detected, image, fwhm=4.0, aper=12.0
    )

    # Match and compute errors
    from scipy.spatial import cKDTree
    tree = cKDTree(np.c_[true_cat['x'], true_cat['y']])
    distances, indices = tree.query(np.c_[measured['x'], measured['y']])

    matched = distances < 3.0
    flux_ratio = measured['flux'][matched] / true_cat['flux'][indices[matched]]
    flux_scatter = np.std(flux_ratio)

    flux_errors.append(flux_scatter)
    print(f"Defocus {defocus} waves: scatter = {flux_scatter*100:.1f}%")

# Plot results
import matplotlib.pyplot as plt
plt.plot(defocus_levels, flux_errors, 'o-')
plt.xlabel('Defocus (waves)')
plt.ylabel('Photometry scatter')
plt.title('Photometry Robustness vs Aberration')
plt.grid(True, alpha=0.3)
plt.show()

Example: Training Real-Bogus Classifier

from stdpipe import simulation
import numpy as np

# Generate diverse training set with various PSF qualities
all_cutouts = []
all_labels = []

for img_idx in range(100):
    # Randomize aberrations
    defocus = np.random.uniform(0.0, 2.0)
    astigmatism = np.random.uniform(0.0, 1.0)
    coma = np.random.uniform(0.0, 0.5)

    # Create aberrated PSF
    psf_model = simulation.create_psf_model(
        fwhm=np.random.uniform(2.5, 5.0),
        psf_type='gaussian',
        defocus=defocus,
        astigmatism_x=astigmatism * np.cos(np.random.uniform(0, 2*np.pi)),
        astigmatism_y=astigmatism * np.sin(np.random.uniform(0, 2*np.pi)),
        coma_x=coma * np.cos(np.random.uniform(0, 2*np.pi)),
        coma_y=coma * np.sin(np.random.uniform(0, 2*np.pi)),
    )

    # Simulate image with aberrated PSF
    result = simulation.simulate_image(
        width=1000,
        height=1000,
        n_stars=100,
        star_psf=psf_model,
        n_cosmic_rays=20,
        n_hot_pixels=30,
        background=1000.0
    )

    # Extract and label cutouts for classifier training
    # ... (detection, matching, cutout extraction code)

# Now you have training data with diverse PSF conditions
# This makes your classifier robust to varying image quality

Technical Notes

Performance: Aberrated PSF creation takes ~50-100ms (one-time cost). Star placement speed is unchanged as it uses the pre-computed PSF.

Accuracy: The Fourier optics approach accurately models:

• Diffraction effects through the aberrated pupil

• Wave interference patterns

• Proper flux conservation

Limitations:

• Does not model chromatic aberrations (wavelength-dependent)

• Assumes monochromatic light at specified wavelength

• Does not model atmospheric scintillation or temporal variations

• High-order aberrations (j>8) not currently supported

References:

• Noll (1976) - “Zernike polynomials and atmospheric turbulence”, J. Opt. Soc. Am.

• Born & Wolf - “Principles of Optics” (aberration theory)

• Mahajan - “Optical Imaging and Aberrations” (validation formulas)

Galaxy Profiles

Creating Custom Galaxy Models

You can create galaxy models with custom Sersic profiles:

from stdpipe.simulation import create_sersic_profile, place_galaxy
import numpy as np

# Create a Sersic profile stamp
profile = create_sersic_profile(
    size=101,
    x0=50.0,
    y0=50.0,
    amplitude=1.0,
    r_eff=10.0,                 # Effective radius (half-light)
    n=1.0,                      # Sersic index
    ellipticity=0.5,            # 0=circular, 1=line
    position_angle=45.0         # Degrees, 0=vertical
)

# Place a galaxy into an image
image = np.zeros((200, 200))
place_galaxy(
    image,
    x0=100.5,
    y0=100.5,
    flux=10000.0,               # Total integrated flux
    r_eff=15.0,
    n=4.0,                      # de Vaucouleurs profile
    ellipticity=0.6,
    position_angle=120.0,
    gain=1.0
)

Common Sersic Profiles

# Exponential disk (spiral galaxy)
cat_spirals = simulation.add_galaxies(
    image, n=10, n_range=(0.8, 1.2), r_eff_range=(10, 20)
)

# Elliptical galaxies (de Vaucouleurs)
cat_ellipticals = simulation.add_galaxies(
    image, n=5, n_range=(3.5, 4.5), r_eff_range=(5, 15),
    ellipticity_range=(0.3, 0.8)
)

# Compact sources (dwarf galaxies)
cat_dwarfs = simulation.add_galaxies(
    image, n=15, n_range=(0.5, 1.0), r_eff_range=(2, 5),
    flux_range=(200, 1000)
)

Real-Bogus Classifier Training

The simulation module is ideal for training real-bogus classifiers because every object is labeled:

from stdpipe import simulation, photometry
import numpy as np

# Generate training set
result = simulation.simulate_image(
    width=2000,
    height=2000,
    n_stars=500,
    star_flux_range=(100, 100000),
    n_galaxies=100,
    n_cosmic_rays=50,
    n_hot_pixels=100,
    n_satellites=5,
    background=1000.0,
    readnoise=10.0,
    verbose=True
)

image = result['image']
true_catalog = result['catalog']

# Detect objects using STDPipe
detected = photometry.get_objects_sep(image, thresh=5.0)

# Match detected objects to truth
from scipy.spatial import cKDTree

# Build KD-tree for truth catalog
truth_coords = np.c_[true_catalog['x'], true_catalog['y']]
tree = cKDTree(truth_coords)

# Match detected objects
detected_coords = np.c_[detected['x'], detected['y']]
distances, indices = tree.query(detected_coords)

# Label detected objects
detected['is_real'] = False
detected['matched_type'] = 'none'

match_radius = 3.0  # pixels
matched = distances < match_radius

detected['is_real'][matched] = true_catalog['is_real'][indices[matched]]
detected['matched_type'][matched] = true_catalog['type'][indices[matched]]

# Now you have labeled training data for ML classifiers
real_sources = detected[detected['is_real'] == True]
artifacts = detected[detected['is_real'] == False]

print(f"Real sources: {len(real_sources)}")
print(f"Artifacts: {len(artifacts)}")
print(f"Unmatched detections: {np.sum(~matched)}")

# Extract features for classifier training
features = np.c_[
    detected['fwhm'],
    detected['a'],              # Semi-major axis
    detected['b'],              # Semi-minor axis
    detected['flux'],
    detected['flags']
]

labels = detected['is_real'].astype(int)

# Train your classifier here...
# from sklearn.ensemble import RandomForestClassifier
# clf = RandomForestClassifier()
# clf.fit(features, labels)

Testing Detection Limits

Simulate images at various flux levels to measure detection completeness:

from stdpipe import simulation, photometry
import numpy as np
import matplotlib.pyplot as plt

# Flux levels to test
flux_levels = np.logspace(2, 4, 20)  # 100 to 10000 ADU
completeness = []

for flux in flux_levels:
    # Simulate field with stars at this flux level
    result = simulation.simulate_image(
        width=1000,
        height=1000,
        n_stars=100,
        star_flux_range=(flux * 0.9, flux * 1.1),
        star_fwhm=3.0,
        background=1000.0,
        readnoise=10.0
    )

    image = result['image']
    true_cat = result['catalog']
    true_stars = true_cat[true_cat['type'] == 'star']

    # Detect objects
    detected = photometry.get_objects_sep(image, thresh=5.0)

    # Match and compute completeness
    from scipy.spatial import cKDTree
    tree = cKDTree(np.c_[true_stars['x'], true_stars['y']])
    distances, _ = tree.query(np.c_[detected['x'], detected['y']])
    n_recovered = np.sum(distances < 3.0)

    completeness.append(n_recovered / len(true_stars))
    print(f"Flux {flux:.0f}: {completeness[-1]*100:.1f}% recovered")

# Plot completeness curve
plt.figure(figsize=(8, 6))
plt.semilogx(flux_levels, completeness, 'o-')
plt.xlabel('Flux (ADU)')
plt.ylabel('Completeness')
plt.title('Detection Completeness vs Flux')
plt.grid(True, alpha=0.3)
plt.ylim(0, 1.1)
plt.show()

Testing Photometry Accuracy

Validate photometry algorithms using simulated sources with known flux:

from stdpipe import simulation, photometry, photometry_measure
import numpy as np

# Simulate isolated stars with known fluxes
result = simulation.simulate_image(
    width=1000,
    height=1000,
    n_stars=50,
    star_flux_range=(5000, 20000),
    star_fwhm=3.5,
    star_psf='gaussian',
    n_galaxies=0,
    n_cosmic_rays=0,
    background=1000.0,
    readnoise=10.0,
    edge=50                     # Keep stars isolated
)

image = result['image']
true_catalog = result['catalog']
true_stars = true_catalog[true_catalog['type'] == 'star']

# Detect and measure
detected = photometry.get_objects_sep(image, thresh=5.0)
measured = photometry_measure.measure_objects(
    detected, image, fwhm=3.5, aper=10.0
)

# Match measured to truth
from scipy.spatial import cKDTree
tree = cKDTree(np.c_[true_stars['x'], true_stars['y']])
distances, indices = tree.query(np.c_[measured['x'], measured['y']])

# Analyze photometry accuracy
matched = distances < 2.0
measured_matched = measured[matched]
true_matched = true_stars[indices[matched]]

# Compare fluxes
flux_ratio = measured_matched['flux'] / true_matched['flux']
flux_bias = np.median(flux_ratio) - 1.0
flux_scatter = np.std(flux_ratio)

print(f"Matched {np.sum(matched)}/{len(true_stars)} stars")
print(f"Photometry bias: {flux_bias*100:.2f}%")
print(f"Photometry scatter: {flux_scatter*100:.2f}%")

# Plot results
import matplotlib.pyplot as plt

plt.figure(figsize=(10, 5))

plt.subplot(1, 2, 1)
plt.plot(true_matched['flux'], measured_matched['flux'], 'o', alpha=0.5)
plt.plot([0, 25000], [0, 25000], 'r--', label='Perfect')
plt.xlabel('True Flux (ADU)')
plt.ylabel('Measured Flux (ADU)')
plt.title('Photometry Accuracy')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(1, 2, 2)
plt.hist(flux_ratio, bins=20, alpha=0.7)
plt.axvline(1.0, color='r', linestyle='--', label='Perfect')
plt.axvline(np.median(flux_ratio), color='g', linestyle='-', label='Median')
plt.xlabel('Measured / True Flux')
plt.ylabel('Count')
plt.title(f'Bias: {flux_bias*100:.2f}%')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Catalog Format

The catalog returned by simulation functions is an Astropy Table with the following columns:

Common Columns (all objects)

• x, y - Image coordinates

• type - Object type: ‘star’, ‘galaxy’, ‘cosmic_ray’, ‘hot_pixel’, ‘satellite_trail’, ‘bad_column’, ‘bad_row’

• is_real - Boolean: True for astronomical sources, False for artifacts

• ra, dec - Sky coordinates (if WCS provided, else NaN)

Star-Specific Columns

• flux - Total flux in ADU

• mag - Instrumental magnitude (-2.5 * log10(flux))

• fwhm - PSF FWHM in pixels

• psf_type - PSF model used (‘gaussian’ or ‘moffat’)

Galaxy-Specific Columns

• flux - Total integrated flux

• r_eff - Effective (half-light) radius in pixels

• sersic_n - Sersic index

• ellipticity - Ellipticity (0=circular, 1=line)

• position_angle - Position angle in degrees

Artifact-Specific Columns

For cosmic rays and satellites:

• length - Track length in pixels

• width - Track width in pixels

• angle - Track angle in degrees

• intensity - Peak intensity

For hot pixels:

• intensity - Pixel intensity

For bad columns:

• position - Column/row index

• bad_type - Type: ‘dead’, ‘hot’, or ‘noisy’

• orientation - ‘vertical’ or ‘horizontal’

API Reference

High-Level Functions

stdpipe.simulation.simulate_image(width, height, n_stars=100, star_flux_range=(100, 10000), star_fwhm=3.0, star_psf='gaussian', star_beta=2.5, n_galaxies=20, galaxy_flux_range=(500, 5000), galaxy_r_eff_range=(3, 15), galaxy_n_range=(0.5, 4.0), galaxy_ellipticity_range=(0.0, 0.7), n_cosmic_rays=5, cosmic_ray_length_range=(5, 50), cosmic_ray_width_range=(1, 3), cosmic_ray_intensity_range=(1000, 10000), cosmic_ray_profile='worm', n_hot_pixels=10, hot_pixel_intensity_range=(500, 5000), n_bad_columns=0, bad_column_type='dead', n_satellites=0, satellite_length_range=(100, 400), satellite_width_range=(2, 8), satellite_intensity_range=(5000, 20000), diffraction_spikes=False, spike_threshold=50000, optical_ghosts=False, ghost_threshold=100000, add_companions=False, companion_fraction=0.2, companion_separation_fwhm=(1.0, 3.0), companion_flux_ratio=(0.5, 1.5), background=1000.0, readnoise=10.0, gain=1.0, edge=10, wcs=None, return_catalog=True, return_masks=False, seed=None, verbose=False)

Simulate a realistic astronomical image with stars, galaxies, and contaminants.

This is the high-level API that creates a complete simulated image in one call.

Parameters:

widthint, optional

Image width in pixels.

heightint, optional

Image height in pixels.

n_starsint, optional

Number of stars to add.

star_flux_rangetuple, optional

(min, max) star flux in ADU.

star_fwhmfloat, optional

FWHM of stellar PSF in pixels.

star_psfstr or dict, optional

PSF type — either a string (‘gaussian’, ‘moffat’) or a PSF model dict (from create_psf_model). Use dict to specify aberrated PSFs.

star_betafloat, optional

Moffat beta parameter (if star_psf=’moffat’).

n_galaxiesint, optional

Number of galaxies to add.

galaxy_flux_rangetuple, optional

(min, max) galaxy flux in ADU.

galaxy_r_eff_rangetuple, optional

(min, max) effective radius in pixels.

galaxy_n_rangetuple, optional

(min, max) Sersic index.

galaxy_ellipticity_rangetuple, optional

(min, max) ellipticity.

n_cosmic_raysint, optional

Number of cosmic ray tracks.

cosmic_ray_length_rangetuple, optional

(min, max) cosmic ray length in pixels.

cosmic_ray_width_rangetuple, optional

(min, max) cosmic ray width in pixels.

cosmic_ray_intensity_rangetuple, optional

(min, max) cosmic ray peak intensity.

cosmic_ray_profilestr, optional

Cosmic ray profile (‘sharp’, ‘tapered’, ‘worm’).

n_hot_pixelsint, optional

Number of hot pixels.

hot_pixel_intensity_rangetuple, optional

(min, max) hot pixel intensity.

n_bad_columnsint, optional

Number of bad columns.

bad_column_typestr, optional

Type of bad columns (‘dead’, ‘hot’, ‘noisy’).

n_satellitesint, optional

Number of satellite trails.

satellite_length_rangetuple, optional

(min, max) satellite trail length in pixels.

satellite_width_rangetuple, optional

(min, max) satellite trail width in pixels.

satellite_intensity_rangetuple, optional

(min, max) satellite trail intensity.

diffraction_spikesbool, optional

Add diffraction spikes to bright stars.

spike_thresholdfloat, optional

Flux threshold for diffraction spikes.

optical_ghostsbool, optional

Add optical ghosts to very bright stars.

ghost_thresholdfloat, optional

Flux threshold for optical ghosts.

add_companionsbool, optional

If True, add close stellar companions to simulate crowded fields.

companion_fractionfloat, optional

Fraction of stars (0-1) that will have companions.

companion_separation_fwhmtuple, optional

(min, max) companion separation in FWHM units.

companion_flux_ratiotuple, optional

(min, max) companion flux relative to primary star.

backgroundfloat, optional

Background level in ADU.

readnoisefloat, optional

Read noise in ADU.

gainfloat, optional

Detector gain in e-/ADU.

edgeint, optional

Minimum distance from image edges for sources.

wcsastropy.wcs.WCS, optional

WCS object for computing sky coordinates.

return_catalogbool, optional

If True, return catalog of all injected sources.

return_masksbool, optional

If True, return separate masks for each artifact type.

seedint, optional

Random seed for reproducibility. If set, calls np.random.seed(seed).

verbosebool, optional

Enable verbose output.

Returns:

dict

Dictionary with ‘image’, ‘catalog’ (if requested), ‘masks’ (if requested), ‘background’, ‘noise’.

Incremental Functions

stdpipe.simulation.add_stars(image, n=100, flux_range=(100, 10000), fwhm=3.0, psf='gaussian', beta=2.5, edge=0, saturation=None, diffraction_spikes=False, spike_threshold=50000, optical_ghosts=False, ghost_threshold=100000, wcs=None)

Add stars to an image with realistic PSF.

Parameters:

imagenumpy.ndarray

Target image (modified in-place).

nint, optional

Number of stars to add.

flux_rangetuple, optional

(min, max) flux range in ADU.

fwhmfloat, optional

FWHM of the PSF in pixels (used if psf is a string).

psfstr or dict, optional

PSF specification — either a string (‘gaussian’, ‘moffat’) or a PSF model dictionary (PSFEx format).

betafloat, optional

Moffat beta parameter (only used if psf=’moffat’).

edgeint, optional

Minimum distance from image edges.

saturationfloat, optional

Saturation level in ADU.

diffraction_spikesbool, optional

If True, add diffraction spikes to bright stars.

spike_thresholdfloat, optional

Flux threshold for adding diffraction spikes.

optical_ghostsbool, optional

If True, add optical ghosts to very bright stars.

ghost_thresholdfloat, optional

Flux threshold for adding optical ghosts.

wcsastropy.wcs.WCS, optional

WCS object for computing RA/Dec.

Returns:

astropy.table.Table

Catalog of added stars.

stdpipe.simulation.add_galaxies(image, n=20, flux_range=(500, 5000), r_eff_range=(3, 15), n_range=(0.5, 4.0), ellipticity_range=(0.0, 0.7), edge=0, wcs=None)

Add galaxies with Sersic profiles to an image.

Parameters:

imagenumpy.ndarray

Target image (modified in-place).

nint, optional

Number of galaxies to add.

flux_rangetuple, optional

(min, max) total flux in ADU.

r_eff_rangetuple, optional

(min, max) effective radius in pixels.

n_rangetuple, optional

(min, max) Sersic index.

ellipticity_rangetuple, optional

(min, max) ellipticity.

edgeint, optional

Minimum distance from image edges.

wcsastropy.wcs.WCS, optional

WCS object for computing RA/Dec.

Returns:

astropy.table.Table

Catalog of added galaxies.

stdpipe.simulation.add_cosmic_rays(image, n_rays=5, length_range=(5, 50), width_range=(1, 3), intensity_range=(1000, 10000), profile='sharp')

Add cosmic ray tracks to an image.

Cosmic rays are placed at random positions with random orientations.

Parameters:

imagenumpy.ndarray

Target image (modified in-place).

n_raysint, optional

Number of cosmic rays to add.

length_rangetuple, optional

(min, max) length in pixels.

width_rangetuple, optional

(min, max) width in pixels.

intensity_rangetuple, optional

(min, max) peak intensity in ADU.

profilestr, optional

Intensity profile (‘sharp’, ‘tapered’, ‘worm’).

Returns:

astropy.table.Table

Catalog of added cosmic rays with columns: x, y, length, width, angle, intensity, type.

stdpipe.simulation.add_hot_pixels(image, n_pixels=10, intensity_range=(500, 5000), clustering=False, cluster_size=3)

Add hot pixels to an image.

Hot pixels are single bright pixels at random locations, optionally clustered to simulate physical detector defects.

Parameters:

imagenumpy.ndarray

Target image (modified in-place).

n_pixelsint, optional

Number of hot pixels to add.

intensity_rangetuple, optional

(min, max) intensity in ADU.

clusteringbool, optional

If True, create clusters of hot pixels.

cluster_sizeint, optional

Number of pixels per cluster (if clustering=True).

Returns:

astropy.table.Table

Catalog of added hot pixels with columns: x, y, intensity, type.

stdpipe.simulation.add_satellite_trails(image, n_trails=1, length_range=(100, 400), width_range=(2, 8), intensity_range=(5000, 20000), profile='linear', tumbling_prob=0.2)

Add satellite trails to an image.

Parameters:

imagenumpy.ndarray

Target image (modified in-place).

n_trailsint, optional

Number of satellite trails to add.

length_rangetuple, optional

(min, max) length in pixels.

width_rangetuple, optional

(min, max) width in pixels.

intensity_rangetuple, optional

(min, max) peak intensity in ADU.

profilestr, optional

Transverse profile (‘linear’, ‘gaussian’).

tumbling_probfloat, optional

Probability of tumbling satellite (intensity variations).

Returns:

astropy.table.Table

Catalog of added trails.

stdpipe.simulation.add_bad_columns(image, n_columns=2, intensity_range=None, bad_type='dead', orientation='vertical')

Add bad columns or rows to an image.

Bad columns can be dead (low/zero values), hot (high values), or noisy (high variance).

Parameters:

imagenumpy.ndarray

Target image (modified in-place).

n_columnsint, optional

Number of bad columns/rows to add.

intensity_rangetuple, optional

(min, max) intensity for ‘hot’ type. Ignored for ‘dead’ and ‘noisy’.

bad_typestr, optional

Type of bad column (‘dead’, ‘hot’, ‘noisy’).

orientationstr, optional

‘vertical’ for columns, ‘horizontal’ for rows.

Returns:

astropy.table.Table

Catalog of bad columns with columns: position, bad_type, orientation, type.

PSF Functions

stdpipe.simulation.create_psf_model(fwhm=3.0, psf_type='gaussian', beta=2.5, size=None, oversampling=2, defocus=0.0, astigmatism_x=0.0, astigmatism_y=0.0, coma_x=0.0, coma_y=0.0, wavelength=5.5e-07, pupil_diameter=1.0)

Create an oversampled PSF model compatible with PSFEx format.

This creates a PSF model that can be used with psf.place_psf_stamp() and other STDPipe routines. Uses oversampling for accurate flux conservation and fast evaluation with sub-pixel positioning.

Optical aberrations can be added via Zernike polynomial coefficients (Noll ordering). When any aberration coefficient is non-zero, a diffraction PSF is computed via Fourier optics and convolved with the atmospheric seeing PSF (Gaussian or Moffat). When all coefficients are zero, the existing analytical path is used unchanged.

Parameters:

fwhmfloat, optional

Full width at half maximum in pixels.

psf_typestr, optional

Type of PSF model (‘gaussian’ or ‘moffat’).

betafloat, optional

Moffat beta parameter (default 2.5, only used for psf_type=’moffat’).

sizeint, optional

Output stamp size in pixels (after downsampling). If None, automatically sized to capture ~99% of PSF flux (>=8*FWHM).

oversamplingint, optional

Oversampling factor (default 2).

defocusfloat, optional

Zernike Z4 defocus coefficient in waves (default 0.0).

astigmatism_xfloat, optional

Zernike Z5 oblique astigmatism coefficient in waves (default 0.0).

astigmatism_yfloat, optional

Zernike Z6 vertical astigmatism coefficient in waves (default 0.0).

coma_xfloat, optional

Zernike Z7 vertical coma coefficient in waves (default 0.0).

coma_yfloat, optional

Zernike Z8 horizontal coma coefficient in waves (default 0.0).

wavelengthfloat, optional

Observation wavelength in meters (default 550e-9, affects diffraction scale).

pupil_diameterfloat, optional

Pupil diameter in meters (default 1.0, affects diffraction scale).

Returns:

dict

Dictionary with PSFEx-compatible structure containing the PSF model.

Galaxy Functions

stdpipe.simulation.create_sersic_profile(size, x0, y0, amplitude=1.0, r_eff=5.0, n=1.0, ellipticity=0.0, position_angle=0.0)

Create a Sersic profile for galaxy simulation.

The Sersic profile is: I(r) = amplitude * exp(-b_n * ((r/r_eff)^(1/n) - 1)) where b_n ≈ 2n - 0.324 for n >= 0.5 (more accurate formula used internally)

Common cases: - n=0.5: Gaussian-like profile - n=1.0: Exponential disk (typical spiral galaxy disk) - n=4.0: de Vaucouleurs profile (typical elliptical galaxy)

Parameters:

sizeint

Size of the profile stamp (must be odd).

x0float

X center position within the stamp.

y0float

Y center position within the stamp.

amplitudefloat, optional

Central surface brightness amplitude.

r_efffloat, optional

Effective radius in pixels (half-light radius).

nfloat, optional

Sersic index (shape parameter).

ellipticityfloat, optional

Ellipticity (0 = circular, 0.5 = moderate, 0.9 = very elongated).

position_anglefloat, optional

Position angle in degrees (0 = vertical, increases CCW).

Returns:

numpy.ndarray

2D array of shape (size, size) with Sersic profile.

stdpipe.simulation.place_galaxy(image, x0, y0, flux, r_eff=5.0, n=1.0, ellipticity=0.0, position_angle=0.0)

Place a galaxy with Sersic profile into an image.

The galaxy stamp is created, scaled to the requested flux, and added to the image. Similar to psf.place_psf_stamp() but for extended sources.

The image is modified in-place.

Parameters:

imagenumpy.ndarray

Target image (2D array), modified in-place.

x0float

X coordinate where to place the galaxy center.

y0float

Y coordinate where to place the galaxy center.

fluxfloat

Total integrated flux in ADU units.

r_efffloat, optional

Effective radius in pixels.

nfloat, optional

Sersic index.

ellipticityfloat, optional

Ellipticity (0 = circular, 1 = line).

position_anglefloat, optional

Position angle in degrees.

Artifact Functions

stdpipe.simulation.create_cosmic_ray(length, width, angle, max_intensity, profile='sharp')

Create a cosmic ray track.

Cosmic rays appear as elongated tracks with sharp edges, oriented randomly.

Parameters:

lengthfloat

Length of the track in pixels.

widthfloat

Width of the track in pixels.

anglefloat

Angle in degrees (0 = horizontal, increases CCW).

max_intensityfloat

Peak intensity in ADU.

profilestr, optional

Intensity profile (‘sharp’, ‘tapered’, ‘worm’).

Returns:

dict

Dictionary with ‘stamp’ (2D array), ‘size’ (stamp dimensions).

stdpipe.simulation.create_satellite_trail(length, width, intensity, angle=None, profile='linear', tumbling=False)

Create a satellite trail.

Satellite trails are long, thin, bright streaks caused by satellites passing through the field.

Parameters:

lengthfloat

Length of the trail in pixels.

widthfloat

Width of the trail in pixels.

intensityfloat

Peak intensity in ADU.

anglefloat, optional

Angle in degrees (0 = horizontal). If None, random.

profilestr, optional

Transverse profile (‘linear’, ‘gaussian’).

tumblingbool, optional

If True, add intensity variations along the trail (tumbling satellite).

Returns:

dict

Dictionary with ‘stamp’ (2D array), ‘angle’ (actual angle used).

stdpipe.simulation.create_diffraction_spikes(size, x0, y0, star_flux, n_spikes=4, spike_length=50, spike_width=2)

Create diffraction spikes for a bright star.

Diffraction spikes are caused by the support structure of telescopes (refractors, SCTs). Typically 4 spikes at 45-degree intervals.

Parameters:

sizeint

Size of the stamp.

x0float

X center (star position).

y0float

Y center (star position).

star_fluxfloat

Flux of the source star (determines spike intensity).

n_spikesint, optional

Number of spikes (typically 4).

spike_lengthfloat, optional

Length of each spike in pixels.

spike_widthfloat, optional

Width of each spike in pixels.

Returns:

numpy.ndarray

2D array with diffraction spikes.

stdpipe.simulation.create_optical_ghost(size, x0, y0, source_flux, ghost_fraction=0.05, offset=(50, 50), blur_sigma=5.0)

Create an optical ghost (reflection artifact).

Optical ghosts are faint, blurred copies of bright sources caused by reflections in the optical system.

Parameters:

sizeint

Minimum size of the stamp (will be expanded if needed to fit ghost).

x0float

X center of original source within stamp.

y0float

Y center of original source within stamp.

source_fluxfloat

Flux of the original source.

ghost_fractionfloat, optional

Fraction of source flux appearing in ghost (typically 0.01-0.1).

offsettuple, optional

(dx, dy) offset of ghost from source in pixels.

blur_sigmafloat, optional

Gaussian blur sigma for the ghost.

Returns:

dict

Dictionary with ‘stamp’ (2D array) and ‘source_pos’ (x, y) giving the source position within the stamp (needed to align stamp onto the image).