Comparison with Source-Extractor (formerly SExtractor)

Source-Extractor (formerly SExtractor) is a popular astronomical image reduction software that runs in the command line. While both photutils and Source-Extractor serve similar purposes, they do so in subtly different ways and have fundamentally different design principles. Nonetheless, the end-result should be extremely similar, regardless of the choice of software, assuming both are correct.

In this notebook, I compare the results of opticam (photutils-based) to Source-Extractor and show that they are equivalent when configured similarly. For simplicity, image calibrations are omitted from this comparison; instead, opticam’s calibration methods are test in a separate notebook.

A test image

For this example, I’ll use an OPTICAM-MX image of V709 Cas:

from pathlib import Path

from astropy.io import fits
import numpy as np

data_path = Path('sextractor_comparison')

data, header = fits.getdata(data_path / 'V709_Cas_image.fits', header=True)
data = np.asarray(data, dtype=np.float64)

print(repr(header))

As we can see, this image was taken in the 2x2 binning mode, and so has a resolution of 1024x1024. Let’s take a look at the image:

from astropy.visualization import simple_norm
from matplotlib import pyplot as plt

fig, ax = plt.subplots()

ax.imshow(data, origin='lower', norm=simple_norm(data, stretch='log'))

plt.show()

Instead of using opticam.Reducer, since we’re only working with a single image, we’ll compute the image background, identify sources, and perform photometry manually (using opticam’s default values).

By default, opticam uses photutils.segmentation.SourceFinder to perform source detection. This requires specifying an npixels parameter that defines the how many pixels need to be above the specified threshold to constitute a source. opticam defaults to a threshold value of 5. Note: this is the threshold above the background RMS rather than the photometric signal-to-noise ratio. For npixels, opticam defaults to a value of 128 for images taken in the 1x1 binning mode, and then reduces this value by the square of the image binning. In this case, npixels would therefore be set to \(128/2^2 = 32\):

import opticam

threshold = 5

finder = opticam.DefaultFinder(npixels=32)

To compute image backgrounds, opticam uses `photutils.background.Background2D <https://photutils.readthedocs.io/en/stable/user_guide/background.html#d-background-and-noise-estimation>`__. By default, Background2D uses a similar algorithm to Source-Extractor to compute the background. Background2D requires specifying a box_size parameter that determines the size of the background mesh used to compute the two-dimensional background. opticam sets this value to 1/32nd of the image width; in this case, box_size would therefore be set to \(1024 / 32 = 32\):

background = opticam.DefaultBackground(box_size=32)

Let’s compute the background of our image using photutils with opticam’s defaults:

bkg = background(data)
data_clean = data - bkg.background

Let’s look at the background subtracted image and the corresponding background mesh to make sure everything looks reasonable:

from matplotlib import pyplot as plt

fig, ax = plt.subplots()

im = ax.imshow(data_clean, origin='lower', norm=simple_norm(data_clean, stretch='log'))

bkg.plot_meshes(ax=ax, marker='.', color='red', alpha=.5, outlines=True, markersize=5,)

fig.colorbar(im)

plt.show()

The mesh looks to be a reasonable size: boxes are larger than individual sources while being small enough to capture variations in the background across the image. Let’s now identify the sources in this background-subtracted image:

cat = finder(data_clean, 5 * bkg.background_rms)

cat.pprint_all()

As we can see, a number of sources have been identified (22, to be specific). Let’s take a look at these sources:

import numpy as np

source_coords = np.asarray([cat['xcentroid'], cat['ycentroid']]).T

semimajor_sigma = np.median(cat['semimajor_sigma'].value)
radius = 5 * semimajor_sigma  # aperture radius for plotting

from matplotlib import pyplot as plt
from matplotlib.patches import Circle

fig, ax = plt.subplots()

im = ax.imshow(data_clean, origin='lower', norm=simple_norm(data_clean, stretch='log'))

for coords in source_coords:
    ax.add_patch(
        Circle(coords, radius=radius, edgecolor='red', facecolor='none')
    )

fig.colorbar(im)

plt.show()

Source detection looks reasonable. Let’s move on to photometry.

For comparison with Source-Extractor, we’ll only look at aperture photometry. To perform aperture photometry, opticam uses the `photutils.aperture <https://photutils.readthedocs.io/en/stable/user_guide/aperture.html>`__ subpackage. This subpackage defines convenience routines for defining apertures and performing aperture photometry. For simplicity, we’ll use a circular aperture with a 5 pixel radius. By default, opticam defines an `EllipticalAperture <https://photutils.readthedocs.io/en/stable/api/photutils.aperture.EllipticalAperture.html#photutils.aperture.EllipticalAperture>`__, so we’ll make it circular by passing the same value to both axes, and setting the orientation of the ellipse to 0 degrees:

photometer = opticam.AperturePhotometer(
    semimajor_axis=5,
    semiminor_axis=5,
    orientation=0.,
    forced=True,  # use forced photometry since we only have one image
    )

results = photometer.compute(
    image=data_clean,
    bias_var=0.,  # ignore bias variance
    dark_var=0.,  # ignore dark variance
    flat_var=0.,  # ignore flat variance
    background_rms=bkg.background_rms,
    cat_coords=source_coords,
    image_coords=None,  # not needed since forced=True
    psf_params=None,  # not needed since aperture size set manually
    read_noise=0.,  # ignore read noise
    )

fluxes = np.asarray(results['flux'])
flux_errs = np.asarray(results['flux_err'])

One important difference between opticam/photutils and Source-Extractor is how apertures handle partially overlapping pixels. By default, photutils (and, by extension, opticam) computes the exact fractional overlap of each pixel with the aperture. In contrast, Source-Extractor divides pixels into 5x5 subpixels, and includes subpixels whose centres fall within the aperture. photutils can be configured to perform the same subpixelling as Source-Extractor, but this is not the default behaviour. As such, opticam will generally result in slightly different flux values when compared to Source-Extractor.

Source-Extractor

To run Source-Extractor, a number of configuration files are required. We therefore need to tweak these configuration files to match the default behaviour of opticam. In particular, we need to edit default.sex to set DETECT_MINAREA to 32, DETECT_THRESH and ANALYSIS_THRESH to 5, PHOT_APERTURES to 10 (since this quantifies the diameter of the aperture), and BACK_SIZE (background mesh size) to 32. The full configuration file used is given below:

with open(data_path / 'opticam_default.sex', 'r') as file:
    print(file.read())

I also edited the default.param configuration file to only save the following columns in the resulting catalogue:

with open(data_path / 'opticam_default.param', 'r') as file:
    print(file.read())

Source-Extractor is a command line tool, so it’s not as easy as photutils to run it in a notebook like this. However, we can use the `os.system() <https://docs.python.org/3/library/os.html#os.system>`__ function to pass the required commands to the system shell:

import os

image_path = data_path.joinpath('V709_Cas_image.fits').resolve()
os.chdir(os.path.join(os.getcwd(), 'sextractor_comparison'))
os.system(f'source-extractor {image_path} -c opticam_default.sex')
os.chdir(os.path.join(os.getcwd(), '..'))  # go back to original directory

We can see that 22 sources were detected, just like photutils. We can also see that the image background was 117.514 with an RMS of 7.8876. Let’s compare these values to those computed by photutils:

print(round(bkg.background_median, 3), round(bkg.background_rms_median, 4))

The background values are extremely similar, as we would hope. Slight differences are likely due to different sigma clipping parameters. Source-Exractor will iteratively clip pixels within a mesh box until all pixels are within :math:`3 sigma of the median value <https://sextractor.readthedocs.io/en/latest/Background.html#background-estimation>`__. By default, photutils (and, by extension, opticam) uses the same sigma clipping threshold, but only performs up to 10 iterations. Since photutils uses an `astropy.stats.SigmaClip <https://docs.astropy.org/en/stable/api/astropy.stats.SigmaClip.html#astropy.stats.SigmaClip>`__ instance to perform sigma clipping, it can be configured to clip until convergence is acheived, as in Source-Extractor, by defining a SigmaClip instance with maxiters=None. However, this is not the default behaviour. As such, opticam will generally result in slightly different background values when compared to Source-Extractor.

Let’s now take a look at the resulting catalogue:

from astropy.table import Table

sextractor_cat = Table.read(data_path / 'test.cat', format="fits")
sextractor_cat.sort('FLUX_APER', reverse=True)

sextractor_cat.pprint()

We can see that we have the coordinates of each detected source as well as the measured flux values and their corresponding errors. Let’s see how the detected sources compare to those of photutils:

from matplotlib.patches import Rectangle

sextractor_coords = np.asarray([sextractor_cat['X_IMAGE'], sextractor_cat['Y_IMAGE']]).T

fig, ax = plt.subplots()

ax.imshow(data_clean, origin='lower', norm=simple_norm(data_clean, stretch='log'))

for coords in source_coords:
    ax.add_patch(
        Circle(coords, radius=radius, edgecolor='red', facecolor='none')
    )

for coords in sextractor_coords:
    ax.add_patch(
        Rectangle(coords - radius, width=2 * radius, height= 2 * radius, edgecolor='white', facecolor='none', zorder=0)
    )

plt.show()

As we can see, source identification appears identical.

Let’s now compare the measured fluxes:

sextractor_fluxes = np.asarray(sextractor_cat['FLUX_APER'])
sextractor_flux_errs = np.asarray(sextractor_cat['FLUXERR_APER'])

fig, ax = plt.subplots()

ax.plot(sextractor_fluxes, sextractor_flux_errs, 'r+', label='Source-Extractor')
ax.plot(fluxes, flux_errs, 'kx', label='opticam')

ax.set_xscale('log')
ax.legend(fontsize='large')

ax.set_ylabel('Flux error [ADU]', fontsize='large')
ax.set_xlabel('Flux [ADU]', fontsize='large')
ax.minorticks_on()
ax.tick_params(which='both', direction='in', right=True, top=True)

plt.show()

Reassuringly, there looks to be very little difference between the two programs. Let’s compare the results more quantitatively.

err_ratio = flux_errs / sextractor_flux_errs
print(f'Average error ratio: {err_ratio.mean()}')
print(f'Max err ratio: {err_ratio.max()}')
print(f'Max err ratio: {err_ratio.min()}')
print()

ratio = fluxes / sextractor_fluxes
print(f'Average ratio: {ratio.mean()}')
print(f'Max ratio: {ratio.max()}')
print(f'Min ratio: {ratio.min()}')

As we can see, opticam (using photutils) performs very similarly to Source-Extractor when both are configured similarly. On average, the difference between opticam/photutils and Source-Extractor was less than 1% in this example, with extremes of \(\sim 2.5\)%. As noted above, this difference is due to a combination of how opticam and Source-Extractor treat partially overlapping pixels differently, and use different sigma clipping parameters.

assert np.isclose(err_ratio.mean(), 1., rtol=0.01, atol=0.)
assert np.isclose(ratio.mean(), 1., rtol=0.01, atol=0.)

That concludes this comparison between opticam and Source-Extractor.