Identifying sources
phoptic uses photutils to find sources in images. This notebook will demonstrate how to define source finders for use with phoptic, as well as explain phoptic’s default behaviour when no source finder is specified.
Test Image
First thing’s first, let’s open an image that contains some sources. For this example, I’ll use one of the images from the Basic Usage tutorial:
[ ]:
from pathlib import Path
import phoptic
out_dir = Path('out')
phoptic.generate_observations(
out_directory=out_dir / 'data',
n_images=3,
circular_aperture=False,
)
[ ]:
from astropy.io import fits
import numpy as np
import os
files = os.listdir(out_dir / 'data')
file = files[0]
with fits.open(out_dir / 'data' / file) as hdul:
print(repr(hdul[0].header))
image = np.array(hdul[0].data)
binning_factor = int(hdul[0].header['BINNING'][0])
[ ]:
from astropy.visualization import simple_norm
from matplotlib import pyplot as plt
fig, ax = plt.subplots(tight_layout=True)
im = ax.imshow(
image,
norm=simple_norm(image, stretch="sqrt"),
origin="lower",
cmap="Greys",
)
ax.set_xlabel("X")
ax.set_ylabel("Y")
plt.show()
Default Finder
phoptic implements two default source finders: Finder (default) and CrowdedFinder (better for crowded fields, but more expensive). Finder does not implement any source deblending, while CrowdedFinder does. Both of these source finders are wrappers for photutils.segmentation.SourceFinder with some added convenience tailored to OPTICAM.
Let’s use phoptic.Finder to identify the sources in the above image:
[ ]:
from phoptic import DefaultFinder
from photutils.segmentation import detect_threshold
# default value for n_pixels is 128 // binning_factor**2
# default value for border_width is 1/16th of the image width
default_finder = DefaultFinder(n_pixels=128 // 4**2, border_width=image.shape[0] // 16)
default_tbl = default_finder(
image,
threshold=detect_threshold(image, n_sigma=5), # 5 sigma detection threshold
)
print(type(default_tbl))
default_tbl.pprint_all()
When calling a DefaultFinder() instance, an astropy.table.QTable instance is returned that is sorted in descending order of source flux. This is how catalogs are defined in phoptic. We can use the table to visualise which sources have been identified and how their fluxes are ranked:
[ ]:
from matplotlib.patches import Circle
fig, ax = plt.subplots(
tight_layout=True,
)
im = ax.imshow(
image,
norm=simple_norm(image, stretch="log"),
origin="lower",
cmap="Greys",
)
for i, row in enumerate(default_tbl):
xc = row['x_centroid']
yc = row['y_centroid']
ax.add_patch(
Circle(
xy=(xc, yc),
radius=5 * row['semimajor_axis'].value,
facecolor='none',
edgecolor='cyan',
lw=1,
)
)
ax.text(
xc * 1.05,
yc * 1.05,
f'{i + 1}',
color='cyan',
)
ax.set_xlabel("X")
ax.set_ylabel("Y")
plt.show()
Custom Source Finders
Defining the Source Finder
Let’s now define a custom source finder. A custom source finder must be a callable that takes two parameters: image and threshold, and returns a QTable instance. image should be an NDArray containing the image data. threshold defines the threshold for source detection in ADU. The simplest way to define a custom source finder is by using the photutils.segmentation.SourceFinder class, which combines source detection and deblending:
[ ]:
%%writefile out/custom_routines.py
from photutils.segmentation import SourceCatalog, SourceFinder
def custom_finder(image, threshold):
finder = SourceFinder(
n_pixels=128 // 4**2, # same as DefaultFinder
deblend=True, # same as DefaultFinder
n_levels=256, # higher than DefaultFinder (more deblending)
contrast=0, # lower than DefaultFinder (more deblending)
progress_bar=False, # disable progress bar (same as DefaultFinder)
)
segment_img = finder(image, threshold)
tbl = SourceCatalog(image, segment_img).to_table()
tbl.sort('segment_flux', reverse=True) # sort in descending order
return tbl
Due to compatibility issues, routines defined inside of IPython notebooks cannot be used by multiprocessing (used by phoptic). To overcome this issue, we have used the %%writefile magic command to write our custom_finder() function to a temporary module called custom_routines.py in our out directory. This is required in Python \(\geq\) 3.14.
To understand how SourceFinder works, and what parameters can be weaked, I refer to the excellent photutils documentation: https://photutils.readthedocs.io/en/stable/user_guide/segmentation.html.
Under-the-hood, phoptic.DefaultFinder also uses photutils.segmentation.SourceFinder, but here we have defined our custom source finder to use different parameter values for better source deblending (though deblending isn’t really needed for our example images). In addition to the different parameter values, our custom finder differs from DefaultFinder in that DefaultFinder automatically omits sources that are close to the edge of an image (where “close” is defined as 1/16th of
the image width, the same size as the default background pixels). For this example, this will not make a difference but it can be important in practise.
SourceFinder returns a SegmentationImage, but phoptic assumes an astropy.table.QTable is returned. Here we have converted the returned SegmentationImage to a QTable by first converting it to a photutils.segmentation.SourceCatalog, and then calling the to_table() method of SourceCatalog.
Let’s now initialise this custom source finder and use it to identify sources in the above image.
[ ]:
from out.custom_routines import custom_finder
custom_tbl = custom_finder(
image,
threshold=detect_threshold(image, n_sigma=5), # 5 sigma detection threshold
)
[ ]:
fig, ax = plt.subplots(
tight_layout=True,
)
im = ax.imshow(
image,
norm=simple_norm(image, stretch="log"),
origin="lower",
cmap="Greys",
)
for i, row in enumerate(custom_tbl):
xc = row['x_centroid']
yc = row['y_centroid']
ax.add_patch(
Circle(
xy=(xc, yc),
radius=5 * row['semimajor_axis'].value,
facecolor='none',
edgecolor='cyan',
lw=1,
)
)
ax.text(
xc * 1.05,
yc * 1.05,
f'{i + 1}',
color='cyan',
)
ax.set_xlabel("X")
ax.set_ylabel("Y")
plt.show()
As we can see, we have once again recovered all six sources. Let’s now see how to use this custom source finder with phoptic.Reducer:
[ ]:
custom_reducer = phoptic.Reducer(
out_directory=out_dir / 'reduced' / 'custom_finder',
data_directory=out_dir / 'data',
finder=custom_finder, # use custom source finder
)
custom_reducer.create_catalogs()
Our custom source finder has been used to successfully identify all six sources. Let’s compare this to using DefaultFinder:
[ ]:
default_reducer = phoptic.Reducer(
out_directory=out_dir / 'reduced' / 'default_finder',
data_directory=out_dir / 'data',
)
default_reducer.create_catalogs()
Unsurprisingly, the default source finder has also identified all six sources. Admittedly, this rather simple example does a poor job of demonstrating why defining a custom source finder may be useful, but hopefully it is clear how custom source finders can be implemented. For more information on defining custom source finders, I again refer to the excellent photutils documentation: https://photutils.readthedocs.io/en/stable/user_guide/segmentation.html. Note that phoptic requires that
source finder routines return an astropy.table.QTable instance when called, while photutils returns a SourceCatalog instance. To convert the SourceCatalog to a QTable, you can use the `to_table() <https://photutils.readthedocs.io/en/stable/api/photutils.segmentation.SourceCatalog.html#photutils.segmentation.SourceCatalog.to_table>`__ method.
That concludes the source finder tutorial for phoptic!