Photometry

Background

Photometry deals with measuring the flux of sources (objects) in astronomical images. The measurement is done usually by subtracting the background and summing all the pixels that belong to the source. The output typically has units of counts or electrons, and this can be converted to some magnitude using a zero point - e.g.,

$Mag = ZP - 2.5 \log_{10}(Counts)$

The zero point can be arbitrary, in this case, we usually refer to the magnitude as instrumental magnitude, or it can be related to some magnitude system (e.g., Vega, AB), which is called calibrated magnitudes. Magnitude calibration is discussed in a different document. In this document, we discuss how to measure flux in 2D astronomical images, and we provide some examples in MATLAB.

Pre-photometry processes

Before we perform photometry, the image needs to go through some basic calibration processes, like bias and dark subtraction, linearization, flat-field correction, and fringe correction. These topics are covered in manuals.topic.BiasFlatCalib Two critical steps required before we perform photometric measurements are to estimate the background level and variance of the background around the source. The background is required in order to subtract the sky contribution from the flux, while the variance is required for any source detection and error estimation. These properties are typically estimated either by: Measuring the (robust) mean and std in an annulus around each source. Measuring the background and variance in small regions of the image and interpolating in between.

The signal-to-noise ratio

Before we discuss the details of how photometry is performed, we need to discuss the signal-to-noise ratio (S/N). The signal-to-noise measures the ratio between signal and noise - in other words, this is the inverse of the relative uncertainty. The relative uncertainty in magnitude units is given by

$\Delta{Mag} = \frac{1.086}{S/N}$

Where the constant 1.086 originates from the logarithmic base of the magnitude system. One can define two types of S/N:

1. S/N for detection $(S/N)_{det}$ - This is the S/N for a detection process. Source detection is a hypothesis testing: in this case, we would like to compare the alternative hypothesis that there is a source in some position q, against the null hypothesis that there is no source at position q (only background). Therefore, in this type of S/N, only the background terms appear in the denominator. I.e., the noise contribution of the source is irrelevant.
2. S/N for a measurement process $(S/N)_{meas}$ - This is the S/N for a measurement of some property of the source (like total flux). Here, the denominator should include all the noise terms that may affect the flux of the source, including the flux of the source itself. It is clear that:

$\Big(\frac{S}{N}\Big){\rm det} \ge \Big(\frac{S}{N}\Big){\rm meas}$

For photometry using CCD or CMOS detectors, the typical terms are:

$\Big(\frac{S}{N}\Big)_{\rm det} = \frac{FAt}{\sqrt{aBAt + aDt + aR^{2}}}$

$\Big(\frac{S}{N}\Big)_{\rm meas} = \frac{FAt}{\sqrt{FAt + aBAt + aDt + aR^{2}}}$

Here $F$ is the flux of the source, $t$ is the exposure time, $B$ is the background level per unit time, $D$ is the same as the background but for the dark current contribution, $R$ is the readout noise, $A$ is the effective collecting area of the telescope, and $a$ is the area of the source (e.g., in pixels). For a Gaussian PSF $a=4\pi\sigma^{2}$, where $\sigma$ is the Gaussian sigma-width. The square root in the denominator is due to the fact that we assume that the noise is Poisson and the photons are not correlated (i.e., the std of a Poisson random variable is given by its square root). When $F>B$ and $F>R^{2}$, we say that the noise is dominated by the source (i.e., the source-noise-dominated regime), when $B>F$ and $B>R^{2}$ we are in the background-noise dominated regime, while when $R^{2}$ is the largest we are in the read-noise-dominated regime.

In reality, the sources have some point spread function (PSF). This means that different pixels containing the source contain different counts, and therefore, one can apply different weights to different pixels. We are usually interested in the weights that maximize the S/N. For Gaussian noise, the S/N for detection, it is possible to show that the optimal weights have the exact shape of the PSF (see Zackay & Ofek 2017a; and Ofek & Zackay 2018 for the Poisson noise case). The optimal S/N (i.e., the maximal possible) for the detection of a Gaussian PSF is given by (Ofek & Ben-Ami 2020):

$\Big(\frac{S}{N}\Big) = \sqrt{\int_{0}^{\infty}{2\pi r dr \frac{F^{2}A_{\rm eff}^{2}t_{\rm E}^{2} }{B A_{\rm eff} t_{\rm E} }\frac{e^{-r^2/(\sigma^{2})} }{4\pi^{2}\sigma^{4}} }} = \frac{F A_{\rm eff} t_{\rm E}}{\sqrt{4\pi\sigma^{2} B A_{\rm eff} t_{\rm E} }} $

Here $\sigma$ is the Gaussian PSF sigma-width, and $A_{\rm eff}$ is the telescope's effective aperture. From this expression, we can also calculate the limiting flux of the telescope for source detection at some (S/N) level:

$F_{\rm lim} = \sqrt{4\pi}\Big(\frac{S}{N}\Big) B^{1/2} \sigma (A_{\rm eff} t_{\rm E})^{-1/2}$

We can also calculate the S/N for source measurments (see derivation in Zackay & Ofek 2017a) of a Gaussian PSF:

$(S/N)^{2}_{\rm meas} = FAt - 2\pi\sigma^{2} Bt \ln{(1+\frac{FA}{2\pi \sigma^{2}B})}$

The optimal way to do photometry

An important question is what is the best way to do photometry? Should we sum the counts in some aperture? or do some weighted summation of the flux? Specifically, we are looking for a method that will maximize the S/N for measurements. Under the simplified assumption that the PSF shape (or the shape of the source) is exactly known, and its position is known, the answer is that one needs to give each pixel a different weight. For a Gaussian noise and in the limit that the background dominates the noise, the optimal weighting function is the shape of the PSF. However, when the source noise dominates the noise, the answer is a top-hat function (similar to what is done in aperture photometry). In some cases, there may be additional complications like the fact that the source position is unknown.

Types of photometry

There are many ways to perform photometric measurements. They mainly differ by the shape of the aperture (i.e., the region over we sum all the pixels), and or the weight we give each pixel (i.e., weighted averaging). Another difference is if we assume we know the exact source position, or we let it be a fitted parameter. Below, we summarize the main types of photometry definitions.

Aperture photometry

In aperture photometry, we sum all the pixels in some radius around the center of the source. This is the simplest form of photometry and in some cases, it is preferred due to its simplicity, or near-optimality for bright sources. A few questions arise: What is the best aperture radius? What should we do about the fact that the image is pixelated, and the sum within a circular aperture is not straightforward to calculate.

Aperture radius

Since the variance of the read noise and background are increasing like the aperture radius squared, while the flux in the source is increasing slower then there is an aperture radius that maximizes the S/N for measurements of aperture photometry. There is usually an optimal aperture radius that will maximize this equation. However, this optimal aperture radius is flux dependent. Typically, a reasonable aperture to use is aperture-radius which equals about 1.25 times the FWHM of the PSF.

Aperture pixelization

The main problem with aperture pixelization is that sources at different sub-pixel positions (which we call here pixel phase) will have a different aperture, and this may affect the photometry. Bickerton & Lupton (2013) present an algorithm to solve this problem using sinc interpolation. The author's preference, also implemented in AstroPack is to apply a sinc interpolation, using sub-pixel fft shift, on the pixelated aperture according to the source sub-pixel position. In this case, the aperture remains pixelated, but we use the same pixelization for all the sources.

PSF photometry

In PSF photometry we fit a PSF shape to a source. The PSF is usually estimated based on bright sources, and the fitting process may have several free parameters, including the source position, and flux, and sometimes the background level may be fitted too. In some cases, for example, when we are interested in measuring faint sources with known positions, we are forcing the position of the source and fitting only the flux -- a process called forced PSF photometry. PSF photometry can be done with several levels of complexity: Single source fitting - Each source is fitted individually. Multi-source fitting - groups of nearby sources are fitted simultaneously in order to overcome some blending issues. Multi-iteration process - in the first iteration bright sources are identified, fitted and subtracted. Next, fainter sources are searched for in the subtracted image, and fitted again. Finally, all the sources found in all iterations may be fitted simultaneously. The multi-iteration scheme is commonly used in dense fields where, it is hard to measure the background, and hence find and measure faint sources.

Convolution-based PSF photometry

Fitting a PSF is roughly equivalent, up to some normalization, to convolving the source image with the PSF. Therefore the cross-correlation of an image with the PSF that is being used for source detection can be used for measuring the flux of the source. This is done by reading the value of the pixel in the filtered image, at the position of the source. This is further discussed and explained in Zackay & Ofek (2017a). An important difference between this step and actual PSF fitting is that filtering, when done for example using fft, return the convolution only at the whole pixel position. This is equivalent to fitting the PSF at the whole pixel position, rather than the actual source position. Therefore, this approach may lead to small bias in the photometry.

Other shapes photometry

Aperture photometry has many variants of non-circular apertures or apertures whose radii are selected according to some algorithm. This may include, isophotal photometry, elliptical apertures, Petrosian magnitude, and more.

Notes on error estimation

Since optical detectors count photons, and many times we assume that these photons are not correlated, the relevant statistic is Poisson Statistics. Since in most cases astronomical images contain many counts per pixel, using the Gaussian statistics is well justified. For example, for source detection, it is justified to use the Gaussian approximation when the background rate is above about 1 count (see Ofek & Zackay 2018). In this case, the variance is equal to the expectancy-value. There may be several complications:

1. In ground-based observations the photons directionality from a single source (i.e., the exact position on the detectors) are correlated. This is due to the Earth's turbulence atmosphere that distorts the wavefront that reaches the telescope (i.e., scintillation). This means that astrometric positional uncertainty can not be estimated using the naive approach of using the square root of the number of photons (see e.g., Ofek 2019).
2. In some cases we neglect contributions from additional sources. Examples include the astrometric uncertainty on the source position that affects the position of the aperture. Since the positional uncertainty is larger (see point 1) than the naive Possion noise estimate, this may have some dramatic effects in some cases (e.g., image subtraction; see Zackay, Ofek & Gal-Yam 2016).
3. In case the image is not Nyquist sampled this may affect the photometric errors (see the section on sampling).

Critical Sampling of images

The Nyquist-Shannon sampling theorem states that a sufficient condition for a sample rate that permits a discrete sequence of samples to capture all the information from a continuous-time signal which is band-limited. By band limited we mean that the Fourier transform of the signal is zero above some frequency - i.e., there is no information above some frequency. If the sampling rate is twice the band-limit rate then we say that we have a critical or Nyquist sampling. This means that a critical sampled sequence can be fully reconstructed (or exactly interpolated). With the aid of additional information, one can fully reconstruct some non-critically sampled sequences. In astronomical images, there is a natural band limit, which is the diffraction limit. Specifically, telescopes do not provide spatial information above the frequency of the diffraction limit. The diffraction limit is . How does this affect photometry? In a non-critically sampled image, without additional information, we can not reconstruct the exact shape of the source and therefore any integration operation we will perform will suffer from inaccuracies. Astronomical detectors on ground-based systems typically have pixels that were matched to the seeing and not to the diffraction limit of the telescope. Therefore, it is rare to meet a detector that is critically sampled. However, in practice, due to the seeing effect, at least in long integrations, the high-frequency information content is small and can be neglected (unless extremely high precision is required). For seeing-limited observations, we find that pixel size which is roughly 1/2.3 of the source FWHM (or about 1-sigma width for a Gaussian PSF) provides good enough sampling for many applications (see example in Ofek 2019).

Photometric measurement routines

Several relevant routines for photometry are implemented in AstroPack. We list them according to three categories:

• Function in the imUtil package - These are usually low-level functions that work on matrices.
• Function in the imProc package - These are functions that works on AstroImage, AstroCatalog, AstroPSF, and ImageComponent objects that should be used by the user for image analysis.
• Function in the ds9 class - The ds9 class is designed for using and manipulating ds9 from within MATLAB. It includes several functions to perform on-the-fly photometry on images, and read catalogs associated with images (i.e., the AstroCatalog in the AstroImage class).

Function in the imUtil package

• The imUtil.kernel2 subpackage contain functions to generate 2D PSF functions, like Gaussian, exponential, annulus, etc.
• The imUtil.background subpackage contain functions to estimate the background and variance in images.
• The imUtil.image.moment2 function can calculate the first and second moment of stamps containing sources, and their aperture photometry.
• The imUtil.sources subpackage contains several photometry routines, like aperPhotCube that performs aperture photometry with fft-shifted aperture, psfPhotCube that fits PSF to stamps with sources, findSources that finds sources in an image using matched filtering, and find_measure_sources that find and measure sources.

Function in the imProc package

• The imProc.sources subpackage contains functions for source detection and measurements, including findSources, findMeasureSources, and psfFitPhot.
• The imProc.psf subpackage contains function that generates and analyzes PSF stamps.
• The imProc.background subpackage contains functions for background and variance estimation.
• The imProc.sources.forcedPhot is performing forced photometry on stationary or moving sources.

Function for quick image photometry (the ds9 package)

• ds9.imexam - An interactive tool for inspection of astronomical images including on-the-fly photometry and reading data on sources from an associated AstroCatalog.

Examples

Here are a few examples of image analysis, including source detection and measurements on astronomical images. We assume that the images are bias subtracted and flat-divided.

Aperture photometry

The following example reads images from the disk into an AstroImage object, measures the image background and variance, and stores it in the AstroImage object. The last command in this example, finds the sources using matched filter and measures their properties, including first and second moments, aperture photometry, and convolution-based PSF photometry. The catalog data is stored in the CatData property in the AstroImage object. The CatData property contains an AstroCatalog object with the found sources catalog.

% Select some fits images
List = io.files.filelist('LAST_*.fits');
% Read a FITS image into an AstroImage object
AI = AstroImage(List);
% or use the shortcut:
AI = AstroImage('LAST_*.fits',RegExp',true);
% Estimate the background and variance, in 256x256 size sub images.
% After this step the Back and Var properties of the AstroImage will be populated.
AI = imProc.background.background(AI, 'SubSizeXY',[256 256]);
% find sources, and measure their properties
% After thus step, the CatData property of the AStroImage will be populated with a catalog of sources.
AI = imProc.sources.findMeasureSources(AI);

Display the found sources

You can mark circles around the positions of the found sources on the image using the ds9 tools:

ds9.open
ds9(AI(1))

PSF photometry

After the source detection and aperture photometry is performed, we can add PSF photometry data to the AstroCatalog associated with the AstroImage:

% Construct PSF
% After this step the PSFData property in the Astrocatalog will be
% populated with a PSF
AI = imProc.psf.constructPSF(AI);
% To see the PSF of the first image
AI(1).PSFData.surface
% Estimate the FWHM of the PSF of the fisrt image
AI(1).PSFData.fwhm
% Perform PSF photometry and add it to the AstroCatalog
AI = imProc.sources.psfFitPhot(AI);

Forced Photometry

Forced photometry is important when we would like to perform photometry on faint targets that may not be detected in all the images. The basic routine for forced photometry is imProc.sources.forcedPhot. This routine starts with an AstroImage containing a calibrated image, MASK image and optionally a catalog and PSF, and performs the following steps:

1. Given stationary or moving source coordinates, it looks for an additional GAIA star near the target.
2. For the source + GAIA reference stars, it performs moment estimation and aperture photometry.
3. Estimate the PSF (if not provided)
4. Perform PSF photometry on the targets.
5. Return a MatchedSources object containing selected properties (e.g., magnitudes, RA, Dec) for each source and each epoch.

Comments: In step 2, the 1st moment is calculated and may shift the user-supplied RA/Dec or X/Y. In order to force the position, set 'MomentMaxIter' to 0 (default). In step 4. The default is that up to 10 iterations are performed, and this may shift the final position of the source by +/-0.5 pixels relative to the initial position. You can control this maximum shift using the 'MaxIter' and 'MaxStep' arguments. For moving sources, set the 'Moving' argument to true, and the number of entries in the 'Coo' argument must match the number of epochs (i.e., the number of images in the AstroImage). Another useful script is pipeline.last.forcedPhotAll - this script will go over an entire directory tree, look for images that contain the requested coordinates and perform forced photometry.

An example for the use of the forced photometry routine:

% load the required images into an AstroImage object
AI = AstroImage('*.fits');
% to expedite the process, you can select the images that contains the source position
% for required source position RA, Dec in degrees:
[AI] = isSkyCooInImage(AI, RA, Dec, [], 'deg');
R = imProc.sources.forcedPhot(AI,'Coo',[RA, Dec],'CooUnits','deg');
% by default this routine will perform forced photometry for the target at RA,Dec,
% but also to all the GAIA stars within 300 arcsec from the target.
% This can be sued to photometrically calibrate the flux
% The output is in MatchedSources object.
% The output can be calibrated using lcUtil.zp_external.
% you can use the 'MomentMaxIter', 'SmallStep', 'MaxStep', 'MaxIter' arguments to control by how much
% the exact position of the source will be allowed to change.
% for more information, use:
help imProc.sources.forcedPhot

Forced photometry for moving sources

The forced photometry routine can be applied to moving sources (e.g., asteroids). First, you have to provide a table of coordinates as a function of time from which the target coordinates will be interpolated. For example:

% generate ephemerids for the required solar system object:
[Cat]=celestial.SolarSys.jpl_horizons('ObjectInd','9804','StartJD',[1 1 2023],'StopJD',[1 3 2023],'StepSizeUnits','h');
R = imProc.sources.forcedPhot(AI,'Coo',Cat, 'Moving',true);

Photometric calibration

The conversion of instrumental flux to calibrated flux (with physical units) is called absolute photometric calibration. In order to achieve this goal we usually calibrate our photometry against some sources in which flux in physical units is known. In order to perform photometric calibration we need to find the zero point which converts the magnitude or flux to calibrated magnitude or flux. However, since the calibrated photometry is given in some band (filter), while our flux is measured in a different filter or slightly different filter (e.g., due to differences in the telescope, detector or atmosphere transmission), the zero point alone is not enough, and we need also some color dependent zero point.

Some basic functions that can be used for photometric calibration include:

• imUtil.calib.simplePhotometricZP - Simple photometric ZP solver. Given the instrumental magnitude and calibrated magnitude fit a ZP including color terms, airmass, position, width: $InstMag-calibMag = ZP + a_{1}*color + a_{2}*color^{2} + ... + a_{3}*AM + a_{4}AMcolor + a_{5}*width + a_{6}*X + a_{7}*Y + ...$. Note that the median of each vector is subtracted before the fit.
• imProc.calib.photometricZP - Calculate an absolute photometric calibration to AstroCatalog. Given an AstroCatalog or AstroImage with a catalog, match the sources against a photometric catalog, and calculate the zero point (ZP) of the catalog. If the input is an AstroImage, the photometric ZP will be calculated only if the WCS.Success is true.

Relative photometry

In relative photometry, our goal is to minimize the scatter in the photometry of all the sources.

The following functions are available:

• lcUtil.zp_meddiff - Estimate light curve ZP based on a simple median of differences. In this method, the ZP is calculated relative to a specific image (specified by RefImInd; default is 1). No color term is calculated. The ZP is calculated by selecting the stars that appear in all images with a mean error below some threshold (MagMagErr). The difference in magnitude between successive images is calculated, the median of mag diff is declared as the ZP.

• lsqRelPhot - Perform relative photometry calibration using the linear least square method. This function solves the following linear problem: m_ij = z_i + M_j + alpha*C + ... (see Ofek et al. 2011). By default will perform two iterations, wherein the second iteration, the errors will be taken from the magnitude vs. std diagram and stars with bad measurements will be removed.