SImMER: A Pipeline for Reducing and Analyzing Images of Stars

The Shane Ao infraRed Camera-Spectrograph (ShARCS; Kupke et al. 2012; Gavel et al. 2014) is a near-infrared camera on the Shane 3 m telescope meant to be used in conjunction with the Shane Adaptive Optics system (ShaneAO) at Lick Observatory. When reducing data from this instrument, SImMER uses a bad pixel map that was last updated at the commissioning of the instrument prior to installation in 2014 (Rosalie McGurk, private communication). The camera has a plate scale of 00333 per pixel, with a field of view of 204 × 204.

Artificial bright sources known as "ghost images" may be present in images produced by this instrument due to secondary reflections within the telescope, but they can be identified by their consistent position angle and angular separation across observations. An example of raw Shane data compared to reduced Shane data is depicted in Figure 2. These data are images of K09203794, a star observed by the Kepler spacecraft (Batalha et al. 2010; Borucki et al. 2010; Bryson et al. 2010; Haas et al. 2010; Jenkins et al. 2010; Koch et al. 2010) and observed on the Shane 3 m telescope on 2018 July 23. These observations were obtained on a cloud-free night with seeing of approximately 1.

Zoom In Zoom Out Reset image size

The Palomar High Angular Resolution Observer (PHARO) is an infrared camera (Hayward et al. 2001) that sits behind the P3K adaptive optics system (Dekany et al. 2013) on the Hale 5.1 m telescope at Palomar Observatory. The instrument offers two plate scales: 0025 per pixel and 0040 per pixel.

There is currently not a bad pixel map specific to this instrument available, so we make use of a more general NaN-filtering method that can be applied to other instruments. Benchmarking this approach against the ShARCS bad pixel map on Shane data reveals that in regions on the detector far away from a source, both bad pixel methods produce nearly identical results. However, images reduced with these two methods are slightly different in regions close to sources. These differences are generally on the order of 200 counts per pixel.

Raw PHARO images are separated into four independent quadrants, each quadrant including 512 × 512 pixels; the four parallel signal chains allow for rapid detector readout (Hodapp et al. 1996). When the data are transferred to FITS files, they are written as separate extensions (Hayward et al. 2001). SImMER includes functions that read and flatten this data into single-exposure FITS files, which can subsequently be passed through the rest of the pipeline in the same manner as ShARCS data.

The simplest way to find the center of an image is to identify the positions of its local maxima. In the interest of utilizing a well-tested, computationally efficient, and easily integrated algorithm, we opt for iterative use of the peak_local_max function provided in the scikit-image library. This function first dilates a target image by applying a maximum filter (sliding a box with a user-determined size (default: 100 pixels) over the image and setting all values within the box to the maximum value within the box). In the process, nearby local image maxima are merged. The size of the box determines the separations within which maxima are considered "nearby." This filtered image is then compared to the original image, and the coordinates at which the original image equals the filtered image are set as the true local maxima. A visual demonstration of this process can be found in the scikit-image documentation. ¹⁷

While the runtime of the peak_local_max function (on the order of one millisecond per image) is amenable to application to multiple images, one complication with its direct application is that a threshold must be set by the user for peak detection. This approach is not preferable when reducing data sets containing on the order of a dozen targets. In order to run this procedure automatically with minimal user intervention, we adjust the peak threshold until only two peak coordinates remain. To arrive at this condition, we implement a binary search algorithm to determine the above-mentioned threshold, reducing the search time from a worst-case time-complexity of ${ \mathcal O }(n)$ to a worst-case time-complexity of ${ \mathcal O }({log}(n))$ for varying the threshold. In practice, this approach can result in a speed-up of more than two orders of magnitude; a single registration time decreases from 8.29 s ± 799 ms to 12 ms ± 168 μs. ¹⁸ We find that we need not place additional constraints on minimum peak distance in the peak-detection algorithm. As a caveat, we anticipate that this algorithm may not work well with all PHARO data, which may contain brightening near the edges of frames. ¹⁹

This method is in general approximate, and it is best suited for exploratory data reductions of large data sets. The results may be discrepant from more robust approaches (such as the "saturated" method; Section 2.4.3).

The quick-look method is very susceptible to error induced by observing conditions. On nights with particularly sub-par seeing, for instance, shot noise may strongly influence the derived center of a source. SImMER therefore offers a number of more image registration schemes that are more robust under noisy conditions.

The first of these image registration schemes is the "empirical PSF" procedure. We assume that changes in the PSF are due largely to differential inter-star (e.g., as controlled by airmass and magnitude) and inter-filter performance. Similarly, we assume that temporal variations in the PSF (e.g., due to varying cloud cover or varying AO performance) are negligible over successive exposures. With these assumptions made, we begin with the given position of the primary source on the image as input to SImMER. ²⁰ We next fit the PSF of the source in a single image using a flexible functional form modeled after a two-dimensional Gaussian:

$$\begin{eqnarray}&&F=\displaystyle \frac{1}{2\pi {s}_{x^{\prime} }{s}_{y^{\prime} }}{e}^{(-x{{\prime} }^{2}{/2s}_{x^{\prime} }^{2}-y{{\prime} }^{2}{/2s}_{y^{\prime} }^{2})},\end{eqnarray} \tag{ 1 }$$

where F is the modeled flux of the source, ${s}_{x^{\prime} }$ controls the standard deviation in the $x^{\prime}$ axis, ${s}_{x^{\prime} }$ controls the standard deviation in the $y^{\prime}$ axis, and $x^{\prime}$ and $y^{\prime}$ are rotated by an angle ϕ from the x and y axes, respectively. $x^{\prime}$ and $y^{\prime}$ are also shifted from the origin by a tunable distance in each axis.

We jointly fit the PSF of the central source across successive exposures as follows. Our fitting procedure, using the emcee package (Foreman-Mackey et al. 2013), produces posteriors distributions for each PSF parameter for each exposure. To assess the joint probability distributions of the shared PSF parameters (i.e., ${s}_{x^{\prime} }$, ${s}_{x^{\prime} }$, and θ) under the assumption of independent measurement, we multiply the posterior distributions of each parameter across successive exposures (e.g., Brogi et al. 2017). True simultaneous fitting would likely require the usage of, e.g., nested sampling (Skilling 2006), which can before better than standard Markov Chain Monte Carlo algorithms for strongly multimodal and/or degenerate parameter spaces (e.g., Feroz & Skilling 2013).

With the PSFs and their centers fit, the source centers of each image are more precisely found with a DAOFind (Stetson 1987) call tuned to the fit PSF. Registration proceeds by stacking images according to their centers. This approach can cleanly accommodate rotations between successive exposures (in the ϕ parameter).

An observer may need to reduce saturated data for a variety of reasons. Some science goals necessitate bright star observations—e.g., a search for very faint stellar companions that are fairly distant (e.g., 5to 10) from a bright (e.g., R = 7 mag) primary in a broadband filter (e.g., J band). For instance, the science case in Hirsch et al. (2021) required such a saturation of bright primaries in order to maximize the observed brightness contrasts. In other cases, an exposure might have inadvertently saturated due to an incorrectly estimated exposure time or variable cloud cover.

For saturated targets, the aforementioned peak-finding algorithm may not operate as intended, even for Shane data lacking fringe brightening. In particular, the saturated region may be large enough that multiple spots in a star's point-spread function (PSF) would match with the dilated image, were the dilation size smaller than the stellar PSF. Even if the dilation size were larger than the stellar PSF, the center of the saturated star's PSF would not necessarily be chosen as the "peak."

With these constraints in mind, we adapt a registration procedure with sub-pixel accuracy detailed in Morzinski et al. (2015) for saturated stars. First, we choose a naive center of the star: the center of the raw image. We then incrementally rotate the image around this center, with rotation angles θ spanning 0° to 360°. Finally, we subtract these rotated images from the original image and sum the residuals from all subtractions. We subsequently change the image's center of rotation prior to repeating the process, as in Figure 3. We compute these residual sums over a search box of centers of rotation, taking the center of rotation corresponding to the lowest summed residuals to be the center of our star.

Mathematically, our approach minimizes the residual function

$$\begin{eqnarray}&&R({x}_{\mathrm{center}},{y}_{\mathrm{center}})=\sum _{\theta =0}^{360}({I}_{\mathrm{rot}}(\theta )-{I}_{0}),\end{eqnarray} \tag{ 2 }$$

varying (x_center, y_center) over a search box.

In the original implementation, Morzinski et al. (2015) rotated the image in 10° increments from 5° to 355°. Reducing the number of rotations to three significantly decreases the runtime of our algorithm; as demonstrated in Figure 4, the runtime increases linearly with the number of rotations, so we increase the speed of our program by a factor of 7. In the interest of symmetry, the rotation angles that we choose are 90°, 180°, and 270°. Additionally, the choice of these rotation angles avoids the need for any interpolation of image data; this approach is advantageous because it both reduces the number of computations performed and removes a source of potential numerical error.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To verify the growth of our runtime, we first run the image registration function 10 times at various numbers of rotations; the standard deviation of these timings is taken to be the error in runtime at that number of rotations. These data are then fit by minimizing the negative log likelihood of our function. Next, we use the emcee (Foreman-Mackey et al. 2013) package to sample the posterior distribution about this likelihood maximum. Using the standard stretch sampler move, we run 32 independent chains over 50,000 steps, achieving convergence by checking the integrated autocorrelation time of the chain (Goodman & Weare 2010) and the Gelman-Rubin statistic (Gelman & Rubin 1992). We then repeat this procedure for a quadratic fit. We find that the Bayes factor favors a linear model of runtime growth to a quadratic one 24:1—which, by the criteria of Kass & Raftery (1995), is "strong" evidence against the quadratic model.

Of the image registration modes provided (quick-look, saturated, and multi-source), we determine the saturated mode to be the most robust, though there is a trade-off with speed; this method can take on the order of 2.5 minutes to execute for a single star, as opposed to the roughly 10 ms required for the quick-look method. We will further develop SImMER to reduce the runtime of this method (see Section 3).

If a user is observing wide stellar binaries of similar apparent magnitudes, measured brightness fluctuations unrelated to astrophysical phenomena (in particular, atmospheric turbulence) may cause one star to appear brighter than the other in one frame and fainter than the other in a successive frame. Left untreated, this effect would result in the previously described algorithms attempting to center successive frames of the same target around different stars, disrupting the image registration process. This problem would be exacerbated in observing strategies that make use of wide dither boxes or if the stars are separated by more than a substantial fraction of the image. For instance, if the target star were centered in the image and the companion near the image's edge, a naive reduction could center some images on the companion star and produce strong edge artifacts in the reduced image.

To rectify this issue, we include in our pipeline a multi-source mode. In this mode, the user selects in each image the star that they would like to designate as the primary. After the user clicks on the rough photocenter of their desired primary, the pipeline performs a search restricted to that region. With this being a decidedly more hands-on approach to image reduction, it should be noted that this mode is only recommended to be used when necessary.