Prepare Data for GalfitS

Image Preparation

Overview

This section documents the imaging preparation pipeline in GalfitS. The process is designed to prepare astronomical images for further analysis through the following steps:

  1. Sky Subtraction – Removing the background sky signal using a polynomial fit.

  2. Image Cut and Mask Generation – Extracting a subimage around the target source and creating a corresponding mask.

  3. Output Writing – Saving the processed cut image, mask image, and sigma image to a FITS file.

Step 1: Sky Subtraction

The first step is to subtract the sky background from your image. This is achieved with the sky_subtraction function from the galfits.images module. Here, a 3rd order polynomial is used to model the background, and sigma-clipped statistics are employed for image normalization and visualization.

Example code:

from galfits import images as IM
from astropy.io import fits
from astropy.stats import sigma_clipped_stats
import matplotlib.pyplot as plt
import gsutils  # assumed to be available for image normalization

# Define file name and load the image
fname = './data/tmass_{0}'.format(band)
img = IM.image(fname+'.fits', unit='cR')
hdu = fits.open(fname+'.fits')
header = hdu[0].header
data_org = hdu[0].data

# Define the output file for the sky-subtracted image
ss_imgname = './data/tmass_{0}_ssimg.fits'.format(band)

# Perform sky subtraction with a 3rd order polynomial
maskplus = img.sky_subtraction(polyfit=True, order=3, filepath=ss_imgname,
                               psfFWHM=3, nsigma=3.0, npixels=5,
                               max_grow_size=100, addgrow=0)

# Calculate sigma-clipped statistics for visualization
sky_mean, sky_median, sky_std = sigma_clipped_stats(data_org, sigma=3.0, maxiters=5)
immin = 1 * sky_std
immax = 5 * sky_std

# Visualize the original image, sky-subtracted image, and the subtracted sky
fig = plt.figure(figsize=(18,6))
ax1 = fig.add_subplot(131)
ax1.imshow(gsutils.normimg(data_org, immin, immax, sky=sky_median, frac=0.6),
           cmap='seismic', origin='lower', vmin=-1, vmax=1, interpolation='nearest')
ax2 = fig.add_subplot(132)
ax2.imshow(gsutils.normimg(data_org - img.data.data, immin, immax, sky=sky_median, frac=0.6),
           cmap='seismic', origin='lower', vmin=-1, vmax=1, interpolation='nearest')
ax3 = fig.add_subplot(133)
ax3.imshow(gsutils.normimg(img.data.data, immin, immax, frac=0.6),
           cmap='seismic', origin='lower', vmin=-1, vmax=1, interpolation='nearest')
plt.show()
Sky Subtraction Example

2MASS H-band image showing the results of sky subtraction.

Step 2: Image Cut and Mask Generation

Following sky subtraction, the next steps are to extract a cutout around the target source and to generate a mask to isolate the source. The functions used are:

  • img_cut: Cuts the image based on the source’s right ascension and declination.

  • generate_cutmask: Creates a source mask by detecting pixels above a threshold, optionally deblending overlapping sources.

Example code:

# Define cutout radius in arcseconds
cutr = 2.5  # unit: arcsec

# Extract the image cutout around the target (ra, dec)
imcut, cp = img.img_cut(ra, dec, cutr)
header = hdu[0].header

# Calculate the sigma image for the cutout
sigimg = img.cal_sigma_image(200, gain=1.)
img.cut_sigma_image = sigimg[cp[2]:cp[3], cp[0]:cp[1]]

# Generate the cut mask image using the specified parameters
img.cut_mask_image = img.generate_cutmask(psfFWHM, data=imcut,
                                           nsigma=2.0, npixels=150,
                                           nlevels=32, contrast=0.001,
                                           deblend=False, growmethod=True,
                                           max_grow_size=15, addgrow=addgrow)

# Write the cut images to a FITS file
outname = './cutimg/tmass_{0}_cut.fits'.format(band)
img.write_imcut(ra, dec, outname, header=header)

# Visualize the cutout with the mask overlay
fig = plt.figure(figsize=(8,8))
sky_mean, sky_median, sky_std = sigma_clipped_stats(imcut, sigma=3.0, maxiters=5)
sky_median = 0.0
immin = 5 * sky_std
immax = imcut.max()
frac = 0.4
plt.imshow(gsutils.normimg(imcut, immin, immax, sky=sky_median, frac=frac),
           cmap='seismic', origin='lower', vmin=-1, vmax=1, interpolation='nearest')
plt.imshow(mask0.astype(float), cmap='Blues', alpha=0.5 * mask0.astype(float),
           origin='lower', interpolation='nearest')
plt.savefig(resupath + '{0}_{1}_cutout.png'.format(tagname, band))
plt.show()
Mask and Image Cut Example

Example of an image cutout with an overlaid mask.

Functions Overview

Several key functions in the galfits.images module are utilized in this pipeline:

  • sky_subtraction(sdmask=None, polyfit=True, order=3, filepath=None, psfFWHM=1.5, nsigma=3.0, npixels=5, max_grow_size=100, addgrow=0, usemedian=False) Performs sky background subtraction using either a polynomial fit (default) or median subtraction. Parameters such as order, psfFWHM, and nsigma control the subtraction process.

  • img_cut(ra, dec, cutsize, cutposition=None, IFS=False, usemaxpix=False, sigma_clipped=True, referpix=None, zoomsize=100, spcorrection=[0.,0.]) Cuts the image around a specified celestial coordinate (ra, dec). Optional parameters allow for fine-tuning the cutout region.

  • generate_cutmask(psfFWHM, data=None, sdmask=None, nsigma=3.0, magnification=3., npixels=5, nlevels=32, contrast=0.001, deblend=False, sky_level=None, source_dia=0.5, max_grow_size=100, growmethod=True, nearbygrow=0, addgrow=0, left_source=True) Generates a mask for the cutout image by detecting sources above a given threshold. It supports various parameters for deblending and controlling the mask growth.

  • write_imcut(ra, dec, outputname, header=None, addpixsc=False, addmodel=False) Saves the processed image cutout along with the mask and sigma images to a FITS file. The output file includes all the relevant image components needed for subsequent modeling.

Conclusion

The image preparation module in GalfitS streamlines the process of preparing astronomical images by automating sky subtraction, image cutouts, and mask generation. This ensures that the data used in subsequent modeling are clean, well-calibrated, and precisely focused on the target sources.

Construct Empirical PSF

Overview

This section describes the process for constructing an empirical Point Spread Function (PSF) using the generate_PSFs function from the galfits.images module. The procedure comprises two main stages:

  1. Selection of Real Point Sources: Identify and clean a catalog of candidate stars from the image. This involves source detection, optional cross-matching with the Gaia catalog, and several filtering steps (flux, geometric, isolation, and central pixel checks).

  2. PSF Fitting: Perform PSF fitting on the cleaned sample of point sources. The fitting process optimizes the empirical PSF model, the source centers within their cutout boxes, and the flux normalizations.

Example Code

# read data
imgname = "./J003515.16+084859.9/ls10/170921_055634_g.fits"
outname = "./J003515.16+084859.9/ls10/psf/170921_055634_g_psf.fits"
with fits.open(imgname) as hdu:
    head = hdu[1].header
    data = hdu[1].data

# creating an image instance and generating epsf
img = IM.image(imgname, hdu=1, unit="cR")
sigmaimg = img.cal_sigma_image(expo, gain=gain, fac1=1., fac2=1.)  # expo & gain need to be set to real values

epsf = img.generate_PSFs(
    # Source detection parameters
    psfFWHM=1.7,
    data=data,
    nsigma=3.0,     # for detect_threshold
    npixels=50,     # for detect_sources

    # Real point source catalog construction
    sigmaimg=sigmaimg,
    crossmatch=True,
    starcord_list=None,   # input list of star coordinates: [[ra, dec], ...]
    fluxthreshould=0.,    # drop ultra-faint sources
    qt=0.1,               # axis ratio threshold for cleaning elongated sources
    sigfwhm=3,
    medianfwhm=11.5,
    size=81,
    isolate=0.3,
    central_tol=2,
    droppsf=[3,6,12,13],
    plot=True,
    figurepath='./',

    # PSF fitting parameters
    notfit=False,
    photutils_method=False,
    oversampling=1,
    a_TV=0,
    a_sky=0,
    learning_rate=2e-5,
    num_steps=80000,
    makestackpsf=False,
    filepaths=outname,
    verbose=True,
)

Selection of Real Point Sources

The generate_PSFs function begins by detecting sources in the image to build an initial, “raw” source catalog. This catalog is generated with a simple signal-to-noise (SNR) cut using $$3sigma$$-clipped statistics, but may include extended objects (galaxies), artifacts (CCD blooming, cosmic rays), and even satellite/plane tracks. Therefore, several parameters are provided to refine the selection:

  1. Cross-matching with Gaia Catalog (Optional): - crossmatch = True

    When enabled, the raw source catalog is cross-matched with the Gaia DR3 catalog (which mainly contains stars). The matching radius is set to $$1.5 times text{psfFWHM}$$.

  2. User-Defined Star Catalog: - starcord_list

    You can provide a list of star coordinates (RA, DEC) to override the raw catalog generated by automatic detection.

  3. Flux Threshold Cut: - fluxthreshould = 0.

    This parameter filters out ultra-faint sources by discarding those with fluxes below a fraction of the mean flux of the detected sources.

  4. Geometric Cuts (Axis Ratio and FWHM): - qt = 0.1

    Ensures that only sources with an axis ratio (minor/major) close to 1 (i.e., round sources) are selected.

    • sigfwhm = 3, medianfwhm = 11.5 These parameters restrict the acceptable range of FWHM values. Sources too large or too small (possibly due to artifacts) are filtered out.

  5. Isolation and Central Pixel Check: - size = 81, isolate = 0.3

    Each source is enclosed in an $$text{size} times text{size}$$ pixel box, and only those with minimal contamination from nearby objects (less than the specified isolate fraction) are retained.

    • central_tol = 2 Discards sources if the central region (a circle with radius $$0.25 times text{size}$$) has too many bad/masked pixels.

  6. Manual Dropping of Specific Sources: - droppsf = [3,6,12,13]

    Allows you to manually remove identified problematic sources based on their index from the source catalog.

Real Point Source Catalog Construction

Fig. 1. Summary of the real point source catalog construction process.

PSF Fitting

After constructing a clean sample of point sources, the function proceeds to fit an empirical PSF. With notfit = False, the fitting process optimizes three main components:

  1. Empirical PSF Model: The resulting PSF is a normalized $$text{size} times text{size}$$ image representing the combined information of all selected point sources.

  2. Center Determination: The precise center of each point source within its box is determined.

  3. Flux Normalization: Each source is assumed to be a scaled version of the PSF, and the algorithm optimizes the scaling factors.

The PSF fitting minimizes a loss function that compares the observed data with the PSF model, incorporating regularization terms (controlled by a_TV and a_sky) to penalize overly sharp features or background variations. The process is iterative, controlled by parameters such as:

  • learning_rate: Determines the convergence speed.

  • num_steps: The total number of fitting iterations.

  • oversampling: Adjusts the spatial resolution of the PSF if a value greater than 1 is provided.

  • photutils_method: Option to use the photutils fitting method; here, the built-in GalfitS method is used (set to False).

The final empirical PSF model, along with the derived source centers and flux normalizations, is saved to the output file specified by filepaths.

Selected Star Candidates

Fig. 2. Selected star candidates used for PSF fitting.

Conclusion

The generate_PSFs function seamlessly integrates source detection, catalog cleaning, and PSF fitting. By adjusting its parameters, you can ensure that only high-quality, isolated point sources contribute to the empirical PSF model, thereby improving the accuracy of subsequent image analysis and modeling tasks.