{ "cells": [ { "cell_type": "markdown", "id": "9e97156f", "metadata": {}, "source": [ "# Comparison with `Source-Extractor` (formerly `SExtractor`)\n", "\n", "`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.\n", "\n", "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](test_error_propagation.ipynb).\n", "\n", "## A test image\n", "\n", "For this example, I'll use an OPTICAM-MX image of V709 Cas:" ] }, { "cell_type": "code", "execution_count": null, "id": "50ed71d1", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:20.860669Z", "iopub.status.busy": "2026-01-23T11:21:20.860545Z", "iopub.status.idle": "2026-01-23T11:21:21.000796Z", "shell.execute_reply": "2026-01-23T11:21:21.000500Z" } }, "outputs": [], "source": [ "from pathlib import Path\n", "\n", "from astropy.io import fits\n", "import numpy as np\n", "\n", "data_path = Path('sextractor_comparison')\n", "\n", "data, header = fits.getdata(data_path / 'V709_Cas_image.fits', header=True)\n", "data = np.asarray(data, dtype=np.float64)\n", "\n", "print(repr(header))" ] }, { "cell_type": "markdown", "id": "b8cd76e2", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": null, "id": "480b10ee", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:21.001914Z", "iopub.status.busy": "2026-01-23T11:21:21.001746Z", "iopub.status.idle": "2026-01-23T11:21:21.480897Z", "shell.execute_reply": "2026-01-23T11:21:21.480608Z" } }, "outputs": [], "source": [ "from astropy.visualization import simple_norm\n", "from matplotlib import pyplot as plt\n", "\n", "fig, ax = plt.subplots()\n", "\n", "ax.imshow(data, origin='lower', norm=simple_norm(data, stretch='log'))\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "0024833f", "metadata": {}, "source": [ "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).\n", "\n", "By default, `opticam` uses [photutils.segmentation.SourceFinder](https://photutils.readthedocs.io/en/stable/user_guide/segmentation.html#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$:" ] }, { "cell_type": "code", "execution_count": null, "id": "9cc29137", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:21.482094Z", "iopub.status.busy": "2026-01-23T11:21:21.481851Z", "iopub.status.idle": "2026-01-23T11:21:23.198571Z", "shell.execute_reply": "2026-01-23T11:21:23.198331Z" } }, "outputs": [], "source": [ "import opticam\n", "\n", "threshold = 5\n", "\n", "finder = opticam.DefaultFinder(npixels=32)" ] }, { "cell_type": "markdown", "id": "7ce08446", "metadata": {}, "source": [ "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$:" ] }, { "cell_type": "code", "execution_count": null, "id": "f0501fa2", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:23.199755Z", "iopub.status.busy": "2026-01-23T11:21:23.199579Z", "iopub.status.idle": "2026-01-23T11:21:23.201292Z", "shell.execute_reply": "2026-01-23T11:21:23.201069Z" } }, "outputs": [], "source": [ "background = opticam.DefaultBackground(box_size=32)" ] }, { "cell_type": "markdown", "id": "b37dd2e5", "metadata": {}, "source": [ "Let's compute the background of our image using `photutils` with `opticam`'s defaults:" ] }, { "cell_type": "code", "execution_count": null, "id": "e51a5b48", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:23.201966Z", "iopub.status.busy": "2026-01-23T11:21:23.201876Z", "iopub.status.idle": "2026-01-23T11:21:23.310615Z", "shell.execute_reply": "2026-01-23T11:21:23.310358Z" } }, "outputs": [], "source": [ "bkg = background(data)\n", "data_clean = data - bkg.background" ] }, { "cell_type": "markdown", "id": "bf6012d2", "metadata": {}, "source": [ "Let's look at the background subtracted image and the corresponding background mesh to make sure everything looks reasonable:" ] }, { "cell_type": "code", "execution_count": null, "id": "3bd52dc3", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:23.311768Z", "iopub.status.busy": "2026-01-23T11:21:23.311682Z", "iopub.status.idle": "2026-01-23T11:21:23.693299Z", "shell.execute_reply": "2026-01-23T11:21:23.693072Z" } }, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "\n", "fig, ax = plt.subplots()\n", "\n", "im = ax.imshow(data_clean, origin='lower', norm=simple_norm(data_clean, stretch='log'))\n", "\n", "bkg.plot_meshes(ax=ax, marker='.', color='red', alpha=.5, outlines=True, markersize=5,)\n", "\n", "fig.colorbar(im)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "f0d9d1b8", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": null, "id": "8c93de36", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:23.694147Z", "iopub.status.busy": "2026-01-23T11:21:23.694057Z", "iopub.status.idle": "2026-01-23T11:21:23.785774Z", "shell.execute_reply": "2026-01-23T11:21:23.785546Z" } }, "outputs": [], "source": [ "cat = finder(data_clean, 5 * bkg.background_rms)\n", "\n", "cat.pprint_all()" ] }, { "cell_type": "markdown", "id": "2ed3881d", "metadata": {}, "source": [ "As we can see, a number of sources have been identified (22, to be specific). Let's take a look at these sources:" ] }, { "cell_type": "code", "execution_count": null, "id": "3fa34d74", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:23.786607Z", "iopub.status.busy": "2026-01-23T11:21:23.786523Z", "iopub.status.idle": "2026-01-23T11:21:23.788280Z", "shell.execute_reply": "2026-01-23T11:21:23.788095Z" } }, "outputs": [], "source": [ "import numpy as np\n", "\n", "source_coords = np.asarray([cat['xcentroid'], cat['ycentroid']]).T\n", "\n", "semimajor_sigma = np.median(cat['semimajor_sigma'].value)\n", "radius = 5 * semimajor_sigma # aperture radius for plotting" ] }, { "cell_type": "code", "execution_count": null, "id": "d77a6d2e", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:23.788935Z", "iopub.status.busy": "2026-01-23T11:21:23.788858Z", "iopub.status.idle": "2026-01-23T11:21:23.937702Z", "shell.execute_reply": "2026-01-23T11:21:23.937428Z" } }, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "from matplotlib.patches import Circle\n", "\n", "fig, ax = plt.subplots()\n", "\n", "im = ax.imshow(data_clean, origin='lower', norm=simple_norm(data_clean, stretch='log'))\n", "\n", "for coords in source_coords:\n", " ax.add_patch(\n", " Circle(coords, radius=radius, edgecolor='red', facecolor='none')\n", " )\n", "\n", "fig.colorbar(im)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "5f8a4254", "metadata": {}, "source": [ "Source detection looks reasonable. Let's move on to photometry.\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": null, "id": "d2eead67", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:23.938669Z", "iopub.status.busy": "2026-01-23T11:21:23.938574Z", "iopub.status.idle": "2026-01-23T11:21:24.080818Z", "shell.execute_reply": "2026-01-23T11:21:24.080547Z" } }, "outputs": [], "source": [ "photometer = opticam.AperturePhotometer(\n", " semimajor_axis=5,\n", " semiminor_axis=5,\n", " orientation=0.,\n", " forced=True, # use forced photometry since we only have one image\n", " )\n", "\n", "results = photometer.compute(\n", " image=data_clean,\n", " bias_var=0., # ignore bias variance\n", " dark_var=0., # ignore dark variance\n", " flat_var=0., # ignore flat variance\n", " background_rms=bkg.background_rms,\n", " cat_coords=source_coords,\n", " image_coords=None, # not needed since forced=True\n", " psf_params=None, # not needed since aperture size set manually\n", " read_noise=0., # ignore read noise\n", " )\n", "\n", "fluxes = np.asarray(results['flux'])\n", "flux_errs = np.asarray(results['flux_err'])" ] }, { "cell_type": "markdown", "id": "7bc8d828", "metadata": {}, "source": [ "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](https://sextractor.readthedocs.io/en/latest/Photom.html#fixed-aperture-flux-flux-aper), 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`.\n", "\n", "## `Source-Extractor`\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": null, "id": "1638447b", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:24.082006Z", "iopub.status.busy": "2026-01-23T11:21:24.081913Z", "iopub.status.idle": "2026-01-23T11:21:24.083719Z", "shell.execute_reply": "2026-01-23T11:21:24.083542Z" } }, "outputs": [], "source": [ "with open(data_path / 'opticam_default.sex', 'r') as file:\n", " print(file.read())" ] }, { "cell_type": "markdown", "id": "369fd7e7", "metadata": {}, "source": [ "I also edited the `default.param` configuration file to only save the following columns in the resulting catalogue:" ] }, { "cell_type": "code", "execution_count": null, "id": "ac900850", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:24.084415Z", "iopub.status.busy": "2026-01-23T11:21:24.084336Z", "iopub.status.idle": "2026-01-23T11:21:24.086009Z", "shell.execute_reply": "2026-01-23T11:21:24.085836Z" } }, "outputs": [], "source": [ "with open(data_path / 'opticam_default.param', 'r') as file:\n", " print(file.read())" ] }, { "cell_type": "markdown", "id": "61233822", "metadata": {}, "source": [ "`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: " ] }, { "cell_type": "code", "execution_count": null, "id": "2e4da5ef", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:24.086688Z", "iopub.status.busy": "2026-01-23T11:21:24.086614Z", "iopub.status.idle": "2026-01-23T11:21:24.113919Z", "shell.execute_reply": "2026-01-23T11:21:24.113594Z" } }, "outputs": [], "source": [ "import os\n", "\n", "image_path = data_path.joinpath('V709_Cas_image.fits').resolve()\n", "os.chdir(os.path.join(os.getcwd(), 'sextractor_comparison'))\n", "os.system(f'source-extractor {image_path} -c opticam_default.sex')\n", "os.chdir(os.path.join(os.getcwd(), '..')) # go back to original directory" ] }, { "cell_type": "markdown", "id": "cadda57c", "metadata": {}, "source": [ "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`:" ] }, { "cell_type": "code", "execution_count": null, "id": "5d5b6031", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:24.115127Z", "iopub.status.busy": "2026-01-23T11:21:24.114927Z", "iopub.status.idle": "2026-01-23T11:21:24.117271Z", "shell.execute_reply": "2026-01-23T11:21:24.116982Z" } }, "outputs": [], "source": [ "print(round(bkg.background_median, 3), round(bkg.background_rms_median, 4))" ] }, { "cell_type": "markdown", "id": "b779ff29", "metadata": {}, "source": [ "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 $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](https://photutils.readthedocs.io/en/stable/user_guide/background.html#d-background-and-noise-estimation). 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`.\n", "\n", "Let's now take a look at the resulting catalogue:" ] }, { "cell_type": "code", "execution_count": null, "id": "065e3bec", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:24.118274Z", "iopub.status.busy": "2026-01-23T11:21:24.118160Z", "iopub.status.idle": "2026-01-23T11:21:24.126508Z", "shell.execute_reply": "2026-01-23T11:21:24.126218Z" } }, "outputs": [], "source": [ "from astropy.table import Table\n", "\n", "sextractor_cat = Table.read(data_path / 'test.cat', format=\"fits\")\n", "sextractor_cat.sort('FLUX_APER', reverse=True)\n", "\n", "sextractor_cat.pprint()" ] }, { "cell_type": "markdown", "id": "726973a0", "metadata": {}, "source": [ "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`:" ] }, { "cell_type": "code", "execution_count": null, "id": "9dd003cd", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:24.127443Z", "iopub.status.busy": "2026-01-23T11:21:24.127350Z", "iopub.status.idle": "2026-01-23T11:21:24.255773Z", "shell.execute_reply": "2026-01-23T11:21:24.255443Z" } }, "outputs": [], "source": [ "from matplotlib.patches import Rectangle\n", "\n", "sextractor_coords = np.asarray([sextractor_cat['X_IMAGE'], sextractor_cat['Y_IMAGE']]).T\n", "\n", "fig, ax = plt.subplots()\n", "\n", "ax.imshow(data_clean, origin='lower', norm=simple_norm(data_clean, stretch='log'))\n", "\n", "for coords in source_coords:\n", " ax.add_patch(\n", " Circle(coords, radius=radius, edgecolor='red', facecolor='none')\n", " )\n", "\n", "for coords in sextractor_coords:\n", " ax.add_patch(\n", " Rectangle(coords - radius, width=2 * radius, height= 2 * radius, edgecolor='white', facecolor='none', zorder=0)\n", " )\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "62abe94a", "metadata": {}, "source": [ "As we can see, source identification appears identical.\n", "\n", "Let's now compare the measured fluxes:" ] }, { "cell_type": "code", "execution_count": null, "id": "37090970", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:24.256682Z", "iopub.status.busy": "2026-01-23T11:21:24.256584Z", "iopub.status.idle": "2026-01-23T11:21:24.762304Z", "shell.execute_reply": "2026-01-23T11:21:24.762096Z" } }, "outputs": [], "source": [ "sextractor_fluxes = np.asarray(sextractor_cat['FLUX_APER'])\n", "sextractor_flux_errs = np.asarray(sextractor_cat['FLUXERR_APER'])\n", "\n", "fig, ax = plt.subplots()\n", "\n", "ax.plot(sextractor_fluxes, sextractor_flux_errs, 'r+', label='Source-Extractor')\n", "ax.plot(fluxes, flux_errs, 'kx', label='opticam')\n", "\n", "ax.set_xscale('log')\n", "ax.legend(fontsize='large')\n", "\n", "ax.set_ylabel('Flux error [ADU]', fontsize='large')\n", "ax.set_xlabel('Flux [ADU]', fontsize='large')\n", "ax.minorticks_on()\n", "ax.tick_params(which='both', direction='in', right=True, top=True)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "708f9de2", "metadata": {}, "source": [ "Reassuringly, there looks to be very little difference between the two programs. Let's compare the results more quantitatively." ] }, { "cell_type": "code", "execution_count": null, "id": "6389e778", "metadata": {}, "outputs": [], "source": [ "err_ratio = flux_errs / sextractor_flux_errs\n", "print(f'Average error ratio: {err_ratio.mean()}')\n", "print(f'Max err ratio: {err_ratio.max()}')\n", "print(f'Max err ratio: {err_ratio.min()}')\n", "print()\n", "\n", "ratio = fluxes / sextractor_fluxes\n", "print(f'Average ratio: {ratio.mean()}')\n", "print(f'Max ratio: {ratio.max()}')\n", "print(f'Min ratio: {ratio.min()}')" ] }, { "cell_type": "markdown", "id": "4a21c060", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "id": "fe47adb4", "metadata": { "execution": { "iopub.execute_input": "2026-01-23T11:21:24.763151Z", "iopub.status.busy": "2026-01-23T11:21:24.763063Z", "iopub.status.idle": "2026-01-23T11:21:24.764701Z", "shell.execute_reply": "2026-01-23T11:21:24.764518Z" }, "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "assert np.isclose(err_ratio.mean(), 1., rtol=0.01, atol=0.)\n", "assert np.isclose(ratio.mean(), 1., rtol=0.01, atol=0.)" ] }, { "cell_type": "markdown", "id": "c3f2321b", "metadata": {}, "source": [ "That concludes this comparison between `opticam` and `Source-Extractor`." ] } ], "metadata": { "kernelspec": { "display_name": "opticam6", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.2" } }, "nbformat": 4, "nbformat_minor": 5 }