Zé Vinícius bio photo

Zé Vinícius

O Captain! My Captain!

Email Twitter Github

Second week of the community bonding period is over!

During that week, I primarily took a detailed reading of the classic Stetson’s 1987 paper, which meticulously describes the algorithms implemented in the DAOPHOT software to perform point spread function photometry on crowded fields.

As stated in the previous post, psf photometry of overlapping sources is feasible with DAOPHOT through the loop FIND-GROUP-NSTAR-SUBTRACT-FIND-GROUP… until no more sources are found.

Right now, photutils already has FIND and SUBTRACT routines implemented. So, I started by sketch an interface of the GROUP routine, which looks like this

def daogroup(starlist, crit_separation=None):
    """
    This is an implementation that follows the DAOGROUP algorithm presented
    by Stetson (1987).

    GROUP divides an entire starlist into sets of distinct, self-contained
    groups of mutually overlapping stars.

    GROUP accepts as input a list of stars and their approximate
    brightenesses relative to a model stellar intensity profile for the frame
    and determines which stars are close enough to be capable of
    adversely influencing each others' profile fits.

    Parameters
    ----------
    starlist : `~astropy.Table` or array-like
        List of stars positions.
        If `~astropy.Table`, columns should be named 'x_0' and 'y_0'.
        Additionally, 'flux_0' may also be provided. 
        If array-like, it should be either (x_0, y_0) or (x_0, y_0, flux_0).
        If 'starlist' only contains x_0 and y_0, 'crit_separation' must be
        provided.
    crit_separation : float (optional)
        Distance, in units of pixels, such that any two stars separated by
        less than this distance will be placed in the same group.
        If None, 'flux_0' must be provided in 'starlist'.
    
    Returns
    -------
    group_starlist : list of `~astropy.Table`
        Each `~astropy.Table` in the list corresponds to a group of mutually
        overlapping starts.

    Notes
    -----
    Assuming the psf fwhm to be known, 'crit_separation' may be set to
    k*fwhm, for some positive real k.

    See
    ---
    `~daofind`
    """

    """
    The GROUP algorithm, as described by Stetson, is as follows.

    Take the first star on `starlist`, then the rest of the list is searched
for stars lying within one critical separation of the star in question. If any
stars are found to lie within one critical separation, they are assigned to
its group. Otherwise, the star in question is isolated and may be fitted by a
single model.
    When the end of the star list has been reached, the stars remaining
unassigned are searched for any lying within the critical separation of the
second star in the group; if such are found, they are assigned to that group.
    This process is iterated until no star remains unassigned.

    Some comments on complexity:
        - For a given star, find all stars lying within one critical
        separation takes at most O(n).
        - For a given group, insert stars takes at most O(n).
    """

On the previous post I also commented about the uncertainties on fitted parameters. On this matter, I went ahead and opened this PR #358.

I also played a bit generating simulated data of isolated sources and using psf_photometry to fit them sequentially. The code that I wrote is as follows

import numpy as np
from astropy.table import Table
from astropy.stats import sigma_clipped_stats
from photutils.datasets import make_random_gaussians
from photutils.datasets import make_noise_image
from photutils.datasets import make_gaussian_sources
from photutils import daofind
from photutils import psf
from photutils import CircularAperture
from matplotlib import rcParams
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
rcParams['image.cmap'] = 'viridis'
rcParams['image.aspect'] = 1  # to get images with square pixels
rcParams['figure.figsize'] = (20,10)
rcParams['image.interpolation'] = 'none'

# generate artificial image
ground_truth = make_random_gaussians(30, [500, 1000], [8, 248], [8, 248],
                                     [2, 2], [2, 2], random_state=12345)
shape = (256, 256)
image = (make_gaussian_sources(shape, ground_truth) +
         make_noise_image(shape, type='poisson', mean=5., random_state=12345))

ground_truth.write('input.html')
# estimate background as the median after sigma clipping the sources
_, bkg, std = sigma_clipped_stats(image, sigma=3.0, iters=5)

# find potential sources with daofind
sources = daofind(image - bkg, threshold=5.0*std, fwhm=4.0)
intab = Table(names=['x_0', 'y_0', 'flux_0'], data=[sources['xcentroid'],
              sources['ycentroid'], sources['flux']])
apertures = CircularAperture((sources['xcentroid'], sources['ycentroid']),
                             r=4.)
# perform fitting
fit_info = ['param_cov']
psf_model = psf.IntegratedGaussianPRF(flux=1, sigma=2.0)
fitted_sources = psf.psf_photometry(image - bkg, intab, psf_model,
                                    fitshape=(8,8),
                                    fitter=fitting.LevMarLSQFitter(),
                                    param_uncert=True)
fitted_sources.write('output.html')
residual_image = psf.subtract_psf(image, psf_model, fitted_sources)

# plot results
plt.subplot(1, 2, 1)
plt.imshow(image, origin='lower', interpolation='nearest')
plt.title('Simulated data')
plt.xlabel('x-position (pixel units)')
plt.ylabel('y-position (pixel units)')
apertures.plot(color='red', lw=1.5, alpha=0.5)
plt.subplot(1, 2, 2)
plt.imshow(residual_image, origin='lower', interpolation='nearest')
plt.title('Residual')
plt.show()

print(fitted_sources)

alt text

For this third week I expect to refine the daogroup API and start a similar sketch for the NSTAR routine. Just to give a glimpse, NSTAR basically accepts as input the list of groups returned by GROUP and, for every group, it creates a compound PSF model which will be fittted to all stars within a given group.

Now, to work! :)