Index: trunk/doc/release.2015/ps1.analysis/analysis.tex =================================================================== --- trunk/doc/release.2015/ps1.analysis/analysis.tex (revision 39819) +++ trunk/doc/release.2015/ps1.analysis/analysis.tex (revision 39820) @@ -96,6 +96,4 @@ \section{INTRODUCTION}\label{sec:intro} -\note{more PS1 background} - The Pan-STARRS Image Processing Pipeline is responsible for the basic analysis of images from the Pan-STARRS telescopes Gigapixel Camera. @@ -115,12 +113,13 @@ An additional constraint on the Pan-STARRS system comes from the high -data rate. PS1 produces typically $\sim 700$ GB per night of imaging -data. These images range from high galactic latitudes to the Galactic -Bulge, so large numbers of measurable stars can be expected in much of -the data. The combination of the high precision goals of the -astrometric and photometric measurements and the high data rate (and a -finite computing budget) mean that the process of detecting, -classifying, and measuring the astronomical objects in the image data -stream in a timely fashion are a significant challenge. +data rate. PS1 produces typically $\sim 500$ exposures per night, +corresponding to $\sim 750$ billion pixels of imaging data. The +images range from high galactic latitudes to the Galactic bulge, so +large numbers of measurable stars can be expected in much of the data. +The combination of the high precision goals of the astrometric and +photometric measurements and the high data rate (and a finite +computing budget) mean that the process of detecting, classifying, and +measuring the astronomical objects in the image data stream in a +timely fashion are a significant challenge. In order to achieve these ambitious goals, the object detection, @@ -154,8 +153,8 @@ automated fashion, does it handle 2D variations well? (P. Stetson) -\item Sextractor : pure aperture measurement with rudimentary - object subtraction. pro: fast, widely used, easy to automate. con: - poor object separation in crowded regions, PSF-modeling is only - beta (psfex), what models are available? (E. Bertin) +\item Sextractor : pure aperture measurement with rudimentary object + subtraction. pro: fast, widely used, easy to automate. con: poor + object separation in crowded regions, PSF-modeling was only in beta, + not widely used at the time. (E. Bertin) \item apphot : IRAF-based aperture photometry. pro: widely used. @@ -172,23 +171,50 @@ \end{itemize} -\note{re-phrase this:} The Pan-STARRS IPP team decided that none of -the existing packages met all of their needs, particularly given the -very challenging goals of the project. We decided to redesign the -photometry analysis from scratch, using the lessons learned from the -existing photometry systems. In the process, the object analysis -software would be written using the data analysis C-code library -written for the IPP, \code{psLib}, and the components of the -photometry code would be integrated into the IPP's mid-level astronomy -data analysis toolkit called \code{psModules}. The result is -'PSPhot', which can be used either as a stand-alone C program, or as -callable set of functions. - -\note{discuss the psphot program varients} - -\begin{verbatim} -Other Varients: -* psphotStack -- 5 filter simultaneous fitting -* psphotFullForce -\end{verbatim} +When the IPP development was starting, the existing photometry +packages either did not meet the level of accuracy required or were +required too much human intervention to be considered for the needs of +PS1. In the case of the SDSS Photo tool, the software was judged to +be too tightly integrated to the architecture of SDSS to be easily +re-integrated into the Pan-STARRS pipelin. A new photometry analysis +package was developed using lessons learned from the existing +photometry systems. In the process, the object analysis software was +written using the data analysis C-code library written for the IPP, +\code{psLib}. Components of the photometry code were integrated into +the IPP's mid-level astronomy data analysis toolkit called +\code{psModules}. The result software, 'PSPhot', can be used either +as a stand-alone C program, or as a set of library functions which may +be integrated into other programs + +The main version of PSPhot is a stand-alone program which is run on a +single image or a group of related images representing the data read +from a camera in a single exposure. The images are expected to have +already been detrended so that pixel values are linearly related to +the flux. The gain may be specific by the configuration system, or a +variance image may be supplied. A mask may also be supplied to mark +good, bad, and suspect pixels. Several variants of psphot have also +been used in the PS1 PV3 analysis. + +The version called PSPhotStack accepts a set of images, each +representing the same patch of sky in a different filter, nominally +the full $grizy$ filter set for the analysis of the PS1 PV3 stack +images, though where insufficient data were available in a given +filter, a subset of these filters was processed as a group. As +discussed in detail below, the PSPhotStack analysis includes the +capability of measuring forced PSF photometry in some filter images +based on the position of sources detected in the other filters. It +also include an option to convolve the set of images to a single, +common PSF size across the filters for the purpose of fixed aperture +photometry. + +A second version of PSPhot used in the PV3 analysis is called +PSPhotFullForce. In this version, a set of image all representing the +same pixels are processed together, with the positions of sources to +be analysed loaded from a supplied file. In this version the +analysis, sources are not discovered -- only the supplied sources are +considered. PSF models are determined for each exposure and the +forced PSF photometry is measured for all sources. A subset of +sources may also be used to measure forced galaxy shape parameters. +As described below, a grid of galaxy models are fitted based on the +supplied guess model. \section{PSPhot Design Goals} @@ -224,8 +250,9 @@ level must be reached for images with 250 mas pixels, implying PSPhot must introduce measurement errors less than 1/50th of a - pixel. \note{the choice of F32 parameters places a numerical limit - of 1e-7 on the accuracy of a pixel relative to the size of a chip - (since a single data value is used for X or Y). For the $4800^2$ - GPC chips, this yields a limit of about 0.25 milliarcsecond.} + pixel. The choice of 32 bit floating point data values for the + source centroids places a numerical limit of 1e-7 on the accuracy of + a pixel relative to the size of a chip (since a single data value is + used for X or Y). For the $4800^2$ GPC chips, this yields a limit + of about 0.25 milliarcsecond. \end{itemize} @@ -244,7 +271,7 @@ \item {\bf Flexible non-PSF models} PSPhot must be able to represent - PSF-like objects as well as non-PSF sources. It must be easy to add - new object models as interesting representations of sources are - invented. + PSF-like objects as well as non-PSF sources (e.g., galaxies). It + must be easy to add new object models as interesting representations + of sources are invented. \item {\bf Clean code base} PSPhot should incorporate a high-degree of @@ -256,7 +283,7 @@ provide the user with methods for assessing the different PSF models. -\item {\bf Careful aperture corrections} PSPhot must carefully measure - and correct for the photometric and astrometric trends introduced by - using analytical PSF models. +\item {\bf Careful systematic corrections} PSPhot must carefully + measure and correct for the photometric and astrometric trends + introduced by using analytical PSF models. \item {\bf User Configurable} PSPhot should allow users to change the @@ -276,5 +303,5 @@ \item {\bf Initial object detection} Smooth, find peaks, measure basic - properties + properties. \item {\bf PSF determination} Select PSF candidates, perform model @@ -288,4 +315,7 @@ properties (aperture or PSF) +\item {\bf Extended Source Analysis} Detailed measurements relevant to + galaxies and/or other extended (non-PSF) sources. + \item {\bf Aperture corrections} Measure the curve-of-growth, spatial aperture variations, and background-error corrections. @@ -296,8 +326,8 @@ PSPhot is highly configurable. Users may choose via the configuration -system which of the above analyses are performed. This may be useful -for testing, but may also allow for specialized use cases. For -example, the PSF model may already be available from external -information, in which case the PSF modeling stage can be skipped. +system which of the above analyses are performed. This is useful for +testing, but also allows for specialized use cases. For example, the +PSF model may already be available from external information, in which +case the PSF modeling stage can be skipped. \subsection{Image Preparation} @@ -343,5 +373,37 @@ circumstance, while a pixel in which persistence ghosts have been subtracted might be useful for detection or even analysis of brighter -sources. \note{can I identify which functions respect which sets of masks} +sources. Table~\ref{tab:mask_values} lists the 16 bit values used for +PS1 mask images, along with their description (see \note{Waters et + al. paper} for additional information). + +\begin{table} +\caption{\label{tab:mask_values} PSPhot / GPC1 Mask Image Pixel Values}\vspace{-0.5cm} +\begin{center} +\begin{tabular}{lcl} +\hline +\hline +{\bf Mask Name} & {\bf Mask Value} & {\bf Description} \\ +\hline + DETECTOR & 0x0001 & A detector defect is present. \\ + FLAT & 0x0002 & The flat field model does not calibrate the pixel reliably. \\ + DARK & 0x0004 & The dark model does not calibrate the pixel reliably. \\ + BLANK & 0x0008 & The pixel does not contain valid data. \\ + CTE & 0x0010 & The pixel has poor charge transfer efficiency. \\ + SAT & 0x0020 & The pixel is saturated. \\ + LOW & 0x0040 & The pixel has a lower value than expected. \\ + SUSPECT & 0x0080 & The pixel is suspected of being bad. \\ + BURNTOOL & 0x0080 & The pixel contain an burntool repaired streak. \\ + CR & 0x0100 & A cosmic ray is present. \\ + SPIKE & 0x0200 & A diffraction spike is present. \\ + GHOST & 0x0400 & An optical ghost is present. \\ + STREAK & 0x0800 & A streak is present. \\ + STARCORE & 0x1000 & A bright star core is present. \\ + CONV.BAD & 0x2000 & The pixel is bad after convolution with a bad pixel. \\ + CONV.POOR& 0x4000 & The pixel is poor after convolution with a bad pixel. \\ + MARK & 0x8000 & An internal flag for temporarily marking a pixel. \\ +\hline +\end{tabular} +\end{center} +\end{table} The variance image, if not supplied is constructed by default from the @@ -367,22 +429,36 @@ Since a typical smoothing or warping operation may introduce correlation between 25 - 100 neighboring pixels, the size of such a -covariance image is prohibitive. In practice, however, there are two -extreme cases which generally are relevant. \note{talk about the - covar matrix for a PSF} - -\subsection{Background (Sky) Model} +covariance image is prohibitive. + +Before sources are detected in the image, a model of the background is +subtracted. The image is divided into a grid of background points +with a spacing of 400 pixels. Superpixels of size $800\times 800$ +pixels are used to measure the local background for each background +grid point, thus over-sampling the background spatial variations by a +factor of 2. In the interest of speed, 10,000 randomly selected +{\em unmasked} pixels in these regions are sampled to determine the +background. \note{flesh out the details here}. Bilinear +interpolation is used to generate a full-resolution image from the grid of +background points, and this image is then subtracted from the science +image. The background image and the background standard deviation +image are kept in memory from which the values of \code{SKY} and +\code{SKY\_SIGMA} are calculated for each object in the output catalog. \subsection{Initial Object Detection} + +\subsubsection{Peak Detection} +\label{sec:peaks} The objects are initially detected by finding the location of local peaks in the image. The flux and variance images are smoothed with a -small circularly symmetric kernel using a two-pass 1D Gaussian -(\note{KEYWORD?}). The smoothed flux and variance images are combined -to generate a significance image in signal-to-noise units -\note{including correction for the covariance, if known}. At this -stage, the goal is only to detect the brighter sources, above a user -defined S/N limit (configuration keyword: \code{PEAK\_NSIGMA}). The -detection efficiency for the brighter sources is not strongly -dependent on the form of this smoothing function. +small circularly symmetric kernel using a two-pass 1D Gaussian. The +smoothed flux and variance images are combined to generate a +significance image in signal-to-noise units, including correction for +the covariance, if known. At this stage, the goal is only to detect +the brighter sources, above a user defined S/N limit (configuration +keyword: \code{PEAKS\_NSIGMA\_LIMIT}). A maximum of +\code{PEAKS\_NMAX} are found at this stage. The detection efficiency +for the brighter sources is not strongly dependent on the form of this +smoothing function. The local peaks in the smoothed image are found by first detecting @@ -397,39 +473,4 @@ the maximum $X$ and $Y$ corners of the region. -\subsection{Footprints} - -\note{need to describe the process of generating the source footprints - and then culling the insignificant peaks} - -\subsubsection{Moments and related} - -\note{disucss the Kron mags} - -\note{this section is wrong: we no longer use S/N clipping, but a - Gaussian window function, chosen based on the measured moment} - -Once a collection of peaks have been identified, basic properties of -the objects are measured. First, the local sky flux is measured -within a square annulus with user-defined dimensions -(\code{INNER\_RADIUS} and \code{OUTER\_RADIUS}), using the sample -median. This local background value is then used to calculate the -object first and second moments within a small user-defined aperture -(\code{MOMENT\_RADIUS}). The first-order moments are a good -representation of the object position, while the second-order moments -are a measure of the object shape. The second-order moments are -somewhat sensitive to the size of the aperture and the accuracy of the -background measurement. The moment calculation is only performed -using pixels which exceed a S/N of 1. If, in the process of -calculating the source moments, the S/N limits reject all but \note{3} -or fewer of the source pixels, the peak is identified as being -suspect, and is not used for further analysis. If the measured -centroid coordinates differ from the peak coordinates be a large -amount (\code{MOMENT\_RADIUS}), then the peak is again identified as -being of poor quality and is rejected. In both of these cases, it is -likely that the `peak' was identified in a region of flat flux -distribution or many saturated or edge pixels. - -\subsubsection{Determination of the Peak Coordinates and Errors} - We use the 9 pixels which include the source peak to fit for the position and position errors. We model the peak of the sources as a @@ -448,32 +489,155 @@ of only 0 or 1, we can greatly simplify the chi-square equation to a square matrix equation with the following values: - -%% fix this: -\begin{verbatim} -| 9 0 0 0 6 6 | C_00 | = \sum F_{i,j} -| 0 6 0 0 0 0 | C_10 | = \sum F_{i,j} x -| 0 0 6 0 0 0 | C_01 | = \sum F_{i,j} y -| 0 0 0 6 0 0 | C_11 | = \sum F_{i,j} x y -| 6 0 0 0 6 4 | C_20 | = \sum F_{i,j} x^2 -| 6 0 0 0 4 6 | C_02 | = \sum F_{i,j} y^2 -\end{verbatim} - -The inverse of the 3x3 matrix terms for $C_{00}$, $C_{20}$, and $C_{02}$ is: -\begin{verbatim} -| +5/9 -1/3 -1/3 | -| -1/3 +1/2 0 | -| -1/3 0 +1/2 | -\end{verbatim} - -The location of the peak is determined from the minimum of the +\[ +\left( \begin{array}{cccccc} +9 & 0 & 0 & 0 & 6 & 6 \\ +0 & 6 & 0 & 0 & 0 & 0 \\ +0 & 0 & 6 & 0 & 0 & 0 \\ +0 & 0 & 0 & 6 & 0 & 0 \\ +6 & 0 & 0 & 0 & 6 & 4 \\ +6 & 0 & 0 & 0 & 4 & 6 \\ +\end{array} \right) +\left( \begin{array}{c} +C_{00}\\ +C_{10}\\ +C_{01}\\ +C_{11}\\ +C_{20}\\ +C_{02}\\ +\end{array} \right) += +\left( \begin{array}{c} +\sum F_{i,j} \\ +\sum F_{i,j} x \\ +\sum F_{i,j} y \\ +\sum F_{i,j} x y \\ +\sum F_{i,j} x^2 \\ +\sum F_{i,j} y^2 \\ +\end{array} \right) +\] + +Inverting the 3x3 matrix terms for $C_{00}$, $C_{20}$, and $C_{02}$, +the location of the peak is determined from the minimum of the bi-quadratic function above, and is given by: - \begin{eqnarray} -Det & = & 4 C_{20} C_{02} - C_{11}^2 \\ -x_{min} & = & (C_{11} C_{01} - 2 C_{02} C_{10}) / Det \\ -y_{min} & = & (C_{11} C_{10} - 2 C_{20} C_{01}) / Det \\ +x_{min} & = & (C_{11} C_{01} - 2 C_{02} C_{10}) D^{-1} \\ +y_{min} & = & (C_{11} C_{10} - 2 C_{20} C_{01}) D^{-1} \\ +D & = & 4 C_{20} C_{02} - C_{11}^2 \end{eqnarray} -\note{error on the peak position} +\subsubsection{Footprints} + +The peaks detected in the image may correspond to real sources, but +they may also correspond to noise fluctuations, especially in the +wings of bright stars. PSPhot attempts to identify peaks which may be +formally significant, but are not locally significant. It first +generates a set of ``footprints'', contiguous collections of pixels in +the smoothed significance image above the detection threshold. These +regions are grown by a small amount to avoid errors on rough edges -- +an image of the footprints is convolved with a disk of radius 3 +pixels. Peaks are assigned to the footprints in which they are +contained (note by definition all peaks must be located in a +footprint). + +For any peak which is not the brightest peak in that footprint it is +possible to reach the brightest peak by following the highest valued +pixels between the two peaks. The lowest pixel along this path is the +{\em key col} for this peak (as used in topographic descriptions of a +mountain). If the key col for a given peak is less than +\code{FOOTPRINT\_CULL\_NSIGMA\_DELTA} (4.0) sigmas below the peak of +interest, the peak is considered to be {\em locally insignificant} and +removed from the list of possible detections. In the vicinity of a +saturated star, the rule is somewhat more agressive as the flat-topped +or structured saturated top of a bright star may appear as multiple +peaks with highly significant cols between them. However, this is an +artifact of the proximity to saturation. In this regime, we require +the col to also be a fixed fraction (5\%) of the saturation below the +peak to avoid being marked as locally insignificant. + +\subsubsection{Centroid and higher-order Moments} + +Once a collection of peaks has been identified, a number of basic +properties of the objects related to the first and second moments are +measured. Below, the second moments are used to select candidate +stellar sources to be used in modeling the PSF. + +In order to measure the moments, it is necessary to define an +appropriate aperture in which the moments are measured. We also apply +a ``window function'', down-weighting the pixels by a Gaussian of size +$\sigma_W$ which is chosen to be large compared to the PSF size. The +choice of the window function $\sigma_W$ and the aperture is an +iterative process: for a given value of $\sigma_W$, the PSF stars will +have a measured value of $\sigma$ which is smaller than the true value +due to the window function. \note{generate examples to illustrate + this}. + +To choose the value of $sigma_W$, we try values of (1, 2, 3, 4.5, 6, +9, 12, 18) pixels. For each of these values, we then select candidate +PSF stars based on the distribution of the measured $\sigma_{x,x}, +\sigma_{y,y}$ values. For each test value of $\sigma_w$, determine +the ratio $f = \frac{\sigma_{x,x} + \sigma{y,y}}{2 \sigma_w}$, i.e., +the ratio of the window size to the observed PSF size. We interpolate +to find a value of $\sigma_W$ for which $f$ is expected to be 0.65. +\note{what is the expected ratio of $\sigma_x$ to the true value?}. +We call this value the \code{MOMENTS\_GAUSS\_SIGMA}. We use an +aperture with a radius of \code{PSF\_MOMENTS\_RADIUS} = 4$\times$ +\code{MOMENTS\_GAUSS\_SIGMA} to select the pixels for the measurement. + +Once \code{PSF\_MOMENTS\_SIGMA} has been determined, moments are +measured as defined below. + +\begin{eqnarray} +x_0 & = & \frac{1}{S} \sum_i (f_i - s_i)x_i w_i \\ +y_0 & = & \frac{1}{S} \sum_i (f_i - s_i)y_i w_i \\ +M_{xx} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^2w_i \\ +M_{xy} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)(y_i - y_0)w_i \\ +M_{yy} & = & \frac{1}{S} \sum_i (f_i - s_i)(y_i - y_0)^2w_i \\ +M_{xxx} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^3w_i / r_i \\ +M_{xxy} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^2(y_i - y_0)w_i / r_i \\ +M_{xyy} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)(y_i - y_0)^2w_i / r_i \\ +M_{yyy} & = & \frac{1}{S} \sum_i (f_i - s_i)(y_i - y_0)^3w_i / r_i \\ +M_{xxxx} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^4w_i / r^2_i \\ +M_{xxxy} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^3(y_i - y_0)w_i / r^2_i \\ +M_{xxyy} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^2(y_i - y_0)^2w_i / r^2_i \\ +M_{xyyy} & = & \frac{1}{S} \sum_i (f_i - s_i)(y_i - y_0)(y_i - y_0)^3w_i / r^2_i \\ +M_{yyyy} & = & \frac{1}{S} \sum_i (f_i - s_i)(y_i - y_0)^4w_i / r^2_i +\end{eqnarray} +where $f_i$ is the flux in a pixel; $s_i$ is the local sky value for +that pixel; $w_i$ is the value of the window function for the pixel; +$S = \sum_i (f_i - s_i) w_i$ is the window-weighted sum of the source +flux, used to re-normalize the moments; $r_i$ is the radius of a +pixel, $\sqrt{(x_i - x_0)^2 + (y_i - y_0)^2}$; The sum is performed +over all pixels in the aperture. For the centroid calculation ($x_0, +y_0$), the peak coordinate (see~\ref{sec:peaks}) is used to define the +aperture and the window function; for higher order moments, the +centroid is used to center the window function. + +If the measured centroid coordinates ($x_0, y_0$) differs from the +peak coordinates be a large amount (\code{MOMENT\_RADIUS}), then the +peak is identified as being of poor quality and is rejected. In +both of these cases, it is likely that the `peak' was identified in a +region of flat flux distribution or many saturated or edge pixels. + +In addition to the moments above, a preliminary Kron radius and flux +are also calculated at this stage. In this analysis, the 1st and +half-radial moments are calculated: +\begin{eqnarray} +M_r & = & \frac{1}{S} \sum_i (f_i - s_i)r_i \\ +M_h & = & \frac{1}{S} \sum_i (f_i - s_i)\sqrt{r_i} +\end{eqnarray} +Note that the window function is not applied in the calculation of +these moments. + +The Kron radius is defined the be 2.5$\times$ the first radial moment. +The Kron flux is the sum of (sky-subtracted) pixel fluxes within the +Kron radius. We also calculate the flux in two related annular +apertures: the Kron inner flux is the sum of pixel values for the +annulus $R_1 < r < 2.5 R_1$, while the Kron outer flux is the sum of +pixel values for $2.5 R_1 < r < 4 R_1$. The first radial moment is +limited at the low and high ends by $R_{\rm min} < M_r < R_{\rm max}$ +where $R_{\rm min}$ is the first radial moment of the PSF stars, or +0.75$\times$ \code{MOMENTS\_GAUSS\_SIGMA} if that cannot be +determined. $R_{\rm max}$ is set to \code{PSF\_MOMENTS\_RADIUS}, the +size of the moments aperture. \subsection{PSF Determination} @@ -491,6 +655,7 @@ \begin{eqnarray} f(x,y) & = & I_o exp (-z) + S \\ - R & = & \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y - - y_o)^2}{2\sigma_y^2} + \sigma_{\rm xy}(x - x_o)(y - y_o) + z & = & \frac{x^2}{2\sigma_x^2} + \frac{y^2}{2\sigma_y^2} + \sigma_{\rm xy} x y \\ + x & = & x_{\rm ccd} - x_o \\ + y & = & y_{\rm ccd} - y_o \\ \end{eqnarray} The object model will have a variety of model parameters, in this case @@ -503,18 +668,19 @@ The point-spread-function (PSF) of an image describes the shape of all unresolved objects in the image. In a typical image, the shape of -point sources is not well described by a single functional form; -rather, the shape will vary as a function of position in the image. -The PSF model therefore must describe the parameter variation as a -function of the position of the object on the image. Note that the -object model consists of a certain number of parameters which are -defined by the PSF model, and another set of parameters which are -independent from object to object. For the case of the elliptical -Gaussian model, the PSF parameters would be the shape terms -($\sigma_x, \sigma_y, \sigma_{\rm xy}$) while the independent -parameters would be the centroid, normalization and local sky values -($x_o, y_o, I_o, S$). PSPhot uses a 2-D polynomial to specify the -variation in the PSF parameters as a function of position in the -image. In the case of the elliptical Gaussian, this implies that the -parameters are each a function of the object centroid coordinates: +point sources is not well described by a single function. Instead, +the shape will vary as a function of position in the image. The PSF +model therefore must describe the parameter variation as a function of +the position of the object on the image. Note that the object model +consists of a certain number of parameters which are defined by the +PSF model, and another set of parameters which are independent from +object to object. For the case of the elliptical Gaussian model, the +PSF parameters would be the shape terms ($\sigma_x, \sigma_y, +\sigma_{\rm xy}$) while the independent parameters would be the +centroid, normalization and local sky values ($x_o, y_o, I_o, S$). +PSPhot uses a 2-D polynomial to specify the variation in the PSF +parameters as a function of position in the image \note{or an + interpolated map}. In the case of the elliptical Gaussian, this +implies that the parameters are each a function of the object centroid +coordinates: \begin{eqnarray} \sigma_x & = & f_1(x,y) \\ @@ -579,7 +745,4 @@ and the PSF model parameters. It also defines a specific order for the 4 independent parameters. - -\note{the code may also require that two of the PSF-like parameters -represent the shape in some way}. \subsubsection{PSF Candidate Object Selection} @@ -626,25 +789,18 @@ parameters are far from the minimization values. PSPhot uses the first and second moments to make a good guess for the centroid and -shape parameters for the PSF models. \note{still true? In order to - minimize the impact of close neighbors, the variance values used in - the fit are enhanced by a fraction of the deviation of the - particular pixel value from the model guess.} Any objects which -fail to converge in the fit are flagged as invalid. - -\note{does the variance enhancement introduce too much bias?} - -\note{discuss the convergence criteria, model parameter guesses} +shape parameters for the PSF models. Any objects which fail to +converge in the fit are flagged as invalid. For the resulting collection of object model parameters, the PSF-dependent parameters of the models are all fitted as a function of -position to a 2-D polynomial. The order of this polynomial is (should -be?) a user-defined parameter. The fitting process for these -polynomials is iterative, and rejects the $3-\sigma$ outliers in each -of three passes. This fitting technique results in a robust -measurement of the variation of the PSF model parameters as a function -of position without being excessively biased by individual objects -which fail drastically. Objects whose model parameters are rejected -by this iterative fitting technique are also marked as invalid and -ignored in the later PSF model fitting stages. +position to a 2-D polynomial. The order of this polynomial is a +user-defined parameter. The fitting process for these polynomials is +iterative, and rejects the $3-\sigma$ outliers in each of three +passes. This fitting technique results in a robust measurement of the +variation of the PSF model parameters as a function of position +without being excessively biased by individual objects which fail +drastically. Objects whose model parameters are rejected by this +iterative fitting technique are also marked as invalid and ignored in +the later PSF model fitting stages. All of the PSF-candidate objects are then re-fitted using the PSF @@ -667,5 +823,7 @@ correction is judged to be the best model. -\subsection{Very Bright Stars} +\subsection{Bright Source Analysis} + +\subsubsection{Very Bright Stars} \note{flesh out} @@ -675,22 +833,4 @@ and subtract a radial profile modeled with a spline (?). -\subsection{PSF vs CR vs Extended} - -\subsection{Bright Source Analysis} - -After a PSF model has been determined, PSPhot performs the analysis of -the bright objects in the image. In this stage, all of the objects -with an estimated signal to noise (based on the moments analysis) -greater than a user-set threshold are analysed and subtracted from the -image. An optional successive stage then finds fainter sources and -measures them as well (see Faint Source Analysis, -Section~\ref{faintsources}). In the bright source analysis stage, two -major varients are available. In the primary version, all objects are -examined (in decending order of brightness) and an appropriate models -is determined for each object which is then subtracted; in the -alternate version, the objects are examined (in decending order of -brightness) and the PSF-like objects subtracted first, then the -extended objects are analysed on a second pass. - \subsubsection{Fast Ensemble PSF Fitting} @@ -698,13 +838,16 @@ convenient to subtract off all of the sources, at least as well as possible at this stage. We make the assumption that all sources are -PSF-like. We also assume their position can be taken as the peak of a -2D quadratic function fitted to the peak pixel and its surrounding 8 -pixels. A single linear fit is used to simultaneously measure all -source fluxes. Since the local sky has been subtracted, this -measurement assumes the local sky is zero. - +PSF-like. If the centroid of the source has been determined, we use +this value for its position; otherwise, we use the interpolated +position of the peak. A single linear fit is used to simultaneously +measure all source fluxes. Since the local sky has been subtracted, +this measurement assumes the local sky is zero. We can write a single +$\chi^2$ equation for this image: \[ -\chi^2 = \sum_{\rm pixels} (F_{x,y} - \sum_{\rm sources} A_i PSF[x,y])^2 +\chi^2 = \sum_{\rm pixels} (F_{x,y} - \sum_{\rm sources} A_i P[x_0,y_0])^2 \] +where $F_{x,y}$ is image flux for each pixel, $P[x_0,y_0]$ is the PSF +model realized at the position of source $i$, and $A_i$ is the +normalization for the source. Minimizing this equation with respect to each of the $A_i$ values @@ -712,18 +855,19 @@ \[ M_{i,j} \bar{A_i} = \bar{F_j}\] where $\bar{A_i}$ is the vector of $A_i$ values, the elements of -$M_{i,j}$ consist of the dot product of the unit-flux PSF for source -$i$ and source $i$, and $\bar{F_j}$ is the dot product of the -unit-flux PSF for source $i$ with the pixels corresponding to source -$i$. The dot products are calculated only using pixels within the +$M_{i,j}$ consist of the dot products of the unit-flux PSF for source +$i$ and source $j$, and $\bar{F_j}$ is the dot product of the +unit-flux PSF for source $j$ with the pixels corresponding to source +$j$. The dot products are calculated only using pixels within the source apertures. Since most sources have no overlap with most other sources, this matrix equation results in a very sparse, mostly diagonal square matrix. The dimension is the number of sources, -likely to be 1000s or 10,000s. Such a matrix does not lend itself to -direct inversion. However, an interative solution quickly yields a -result with sufficient accuracy. In the iterative solution, a guess -at the solution is made; the guess is multiplied by the matrix, and -the result compared with the observed vector $\bar{F_j}$. The -difference is used to modify the initial guess. This proces is -repeated several times to achieve a good convergence. +likely to be 1000s or 10,000s. Direct inversion of the matrix would +be computationally very slow. However, an interative solution quickly +yields a result with sufficient accuracy. In the iterative solution, +a guess at the solution $\bar{A}$ is made assuming $M_{i,j}$ is purely +diagonal; the guess is multiplied by $M_{i,j}$, and the result +compared with the observed vector $\bar{F_j}$. The difference is used +to modify the initial guess. This proces is repeated several times to +achieve a good convergence. Once a solution set for $A_i$ is found, all of the objects are @@ -744,5 +888,5 @@ quality of the PSF model as a representation of the object shape. To do this, it calculates the next step of the minimization {\em allowing -the shape parameters to vary}. This step, essentially the + the shape parameters to vary}. This step, essentially the Gauss-Newton minimization distance from the current local minimum, should be very small if the object is well represented by the PSF, but @@ -752,14 +896,13 @@ elliptical Gaussian, these two parameters are $\sigma_x$ and $\sigma_y$. For a generic model, the shape parameters may be defined -differently, but the should always be two parameters which scale the -object size in two dimensions (what about a polar-coordinate form?) -Currently, PSPhot requires the two relevant shape parameters to be the -first two dependent parameters in the list of model parameters (ie, -parameters 4 \& 5). +differently, but there should always be two parameters which scale the +object size in two dimensions. Currently, PSPhot requires the two +relevant shape parameters to be the first two dependent parameters in +the list of model parameters (ie, parameters 4 \& 5). The expected distribution of the variation of the two shape parameters will be a function of the signal-to-noise of the object in question and the value of the shape parameter itself. The expected standard -deviation on the shape parameter is, eg, $\sigma_x / {\rm SN}$. If +deviation on the shape parameter is, eg, $\sigma_x / {\rm S/N}$. If the object is well-represented by the PSF, then the shape parameter values should be close to their minimization value. We can thus ask, @@ -802,7 +945,7 @@ distribution, the remaining flux should be below 1 $\sigma$ significance. In practice the cores of bright stars are poorly -represented and may have larger residual significance. \note{in future -work, we may choose to enhance the variance to minimize detection of -objects in the residuals of brighter objects}. +represented and may have larger residual significance. + +\note{I am not sure the above discussion is still (PV3) true. To be reviewed.} \subsubsection{Blended Sources} @@ -828,6 +971,35 @@ available non-PSF model or models. -\note{better description of the acceptance criteria; the FLT model is - tried before the DBL is accepted or rejected}. +\subsubsection{Source Size Assessment} + +After the PSF model has been fitted to all sources, and the Kron flux +has been measured for all sources, PSPhot uses these two measurements, +along with some additional pixel-level analysis, to determine the size class +of the object. If the object is large compared to a PSF, it is +considered to be {\em extended} and will be +fitted with a galaxy model (or possibly another type of extended +source model in special cases). If the object is small compared to a +PSF, it is considered to be a {\em cosmic ray} and masked. + +Extended sources are identified as those for which the Kron magnitude is +significantly brighter than the PSF magnitude when compared to a PSF +star. The value $dMagKP = m_{\rm Kron} - m_{\rm PSF}$, the difference between the PSF +and Kron magnitudes, is calculated for each object. The median of +$dMagKP$ is calculated for the PSF stars. This median is subtracted +from $dMagKP$ for each star. The result is divided by the quadrature +error of the PSF and Kron magnitudes and called \code{extNsigma}. If +\code{extNsigma} is larger than \code{PSPHOT.EXT.NSIGMA.LIMIT} (3.0), +the object is considered to be extended. + +Cosmic Rays are identified by a combination of the Kron magnitude and +the second-moment width of the object in the minor axis direction. +The second-moment in the minor axis direction is calculated from +$M_{xx}, M_{xy}, M_{yy}$ as follows: +\[ +M_{\rm minor} = \frac{1}{2}(M_{xx} + M_{yy}) - \frac{1}{2}\sqrt{(M_{xx} - M_{yy})^2 + 4 M_{xy}^2} +\] +If $M_{\rm minor} < 1.2$ pixels$^2$ and the instrumental Kron +magnitude is $< -5.5$, then the object is identified as a cosmic ray +and the associated pixels are masked. \subsubsection{Non-PSF Objects} @@ -873,7 +1045,4 @@ \subsection{Faint Sources} -\note{this is not done : should use the ensemble PSF fitting to fit - just the new significant peaks} - After a first pass through the image, in which the brighter sources above a high threshold level have been detected, measured, and @@ -887,57 +1056,9 @@ The objects which are measured in this faint-object stage are clearly low significance detections. A typical threshold for the bright -object analysis would S/N of 5 - 10. The lower limit cutoff for the -faint object analysis would typically be S/N of 2 - 4. In this stage, -PSPhot can perform one of three types of analysis. The difference -between these options is one of speed vs detail. - -In the first option, PSPhot can repeat the analysis described above in -sections XXX and XXX, performing a PSF fit followed by a non-PSF fit -to the objects which are not PSF-like, and subtracting them. The -advantage of this option is that the faint objects are treated -identically to the bright objects, and all potential parameters are -measured, even for marginally extended sources. The disadvantage of -this option is that the low signal-to-noise of the objects in this -stage limits the amount of information which can be extracted from -them. The marginal gain may not be worth the large expense of -processing time. - -In the second option, PSPhot can perform only the PSF model fit to the -remaining peaks, but ignore any further questions of the shape of the -objects. The advantage of this option is that it is substantially -faster than performing the more complex (and less stable) -multi-parameter non-linear fits on all faint objects. On the -downside, less information is learned about the objects. - -Finally, PSPhot can simply ignore the fitting process and instead -glean information about the fainter sources on the basis of the peak -characteristics. In this option, the image is smoothed with the PSF -model, and the peak for each object is measured. The peak flux and -the local peak curvature theoretically give sufficient information to -recover the object flux, the centroid coordinates, and the centroid -errors. The advantage of the stage is speed, especially for the very -faintest levels: if the lower limit is not sufficiently faint, the -time spent in performing the smoothing (3 FFTs) cannot make up for the -time which would have been spent applying the PSF model to the peaks. -The downside of this method is an increased sensitivity to the local -sky model (the local sky must be correctly subtracted) and fewer -constraints on the quality of the detection (no Chi-Square is -measured, for example). - -\note{In the ideal case, if we were only interested in detecting PSFs, -and we had a good model for the PSF, we could optimally find the -sources by smoothing the image and the variance image with the PSF model. -\em write out the description of Nick's optimal PSF finding}. - -PSPhot allows the user to select between these three options for the -analysis of the faint sources. Three separate user-controlled -signal-to-noise ratio limits are defined. One specifies the depth to -which the PSF / non-PSF analysis is performed. A second (which must -be smaller) specifies the depth to which only the PSF is fitted. A -third specifies the depth to which the analysis is performed using on -the peak statistics. If two of these are identical, then certain of -these options are simply skipped. For example, if the peak analysis -threshold is set to the same value as the PSF-only threshold, no peak -analysis is performed. +object analysis would S/N of 5 - 10. \note{PV3 value is 20.0?} The +lower limit cutoff for the faint object analysis would typically be +S/N of 2 - 4. \note{PV3 value is 5.0?} Objects detected in the faint +object stage are fitted with the PSF model using the linear, ensemble +fitting process. \subsection{Aperture Correction Measurement} @@ -990,20 +1111,17 @@ saturation. -\note{this discussion sucks: put in some more details of my point: - amplitude of systematic vs random sky errors} - -How important is this effect? Consider a typical bright object with a -flux of (say) 40,000 counts in an image of background 1000 counts per -pixel, with FWHM of 4 pixels. In principle, the flux of this object -should be measurable with an accuracy of roughly 0.57\% -($\frac{\sqrt{40000 + 1000 \times 12}}{40000}$). However, the -measurement of the sky is limited at some finite level by Poisson -statistics. If we are required to use an aperture of (say) 25 pixels -in radius (eg, 5 arcseconds for an 0.2 arcsec / pixel detector), and -we have an annulus of twice this radius to measure the local sky, then -we will have an error of XXX. - -\note{outline the variation of {\em ApResid} as a function of -magnitude}. +% How important is this effect? Consider a typical bright object with a +% flux of (say) 40,000 counts in an image of background 1000 counts per +% pixel, with FWHM of 4 pixels. In principle, the flux of this object +% should be measurable with an accuracy of roughly 0.57\% +% ($\frac{\sqrt{40000 + 1000 \times 12}}{40000}$). However, the +% measurement of the sky is limited at some finite level by Poisson +% statistics. If we are required to use an aperture of (say) 25 pixels +% in radius (eg, 5 arcseconds for an 0.2 arcsec / pixel detector), and +% we have an annulus of twice this radius to measure the local sky, then +% we will have an error of XXX. +% +% \note{outline the variation of {\em ApResid} as a function of +% magnitude}. PSPhot measures the aperture correction ({\em ApResid}) for every PSF @@ -1325,15 +1443,4 @@ using an FFT-based convolution \note{(examples?)} -Recipe parameters which affect the PSF-convolved galaxy model fitting -process: -\begin{verbatim} -EXT_FIT_NSIGMA_CONV [9] : number of sigma -EXT_FIT_ITER -EXT_FIT_MIN_TOL -EXT_FIT_MAX_TOL -LMM_FIT_CHISQ_CONVERGENCE -LMM_FIT_GAIN_FACTOR_MODE -\end{verbatim} - For the Exponential and DeVaucouleur fits, all parameters are fitted in the non-linear minimization stage. For the Sersic model, we do not @@ -1368,21 +1475,37 @@ \subsection{Convolved Radial Aperture Photometry} +For some science goals, a well-measured color of a galaxy is more +important than an accurate total magnitude. In the case of PS1, the image +quality variations for stacks of different filters presents a serious +challenge for the determination of precise colors. PSPhot determines +a set of PSF-matched radial aperture flux measurements in order to +minimize the impact of the stack image quality variations. + +In PSPhotStack, the stack analysis version of PSPhot, the 5 filter +images are processed together. After the PSF models have been fitted +and a best set of galaxy models have been determined, three sets of +radial apertures are measured. In the first set, the fluxes in the +radial apertures are measured using the raw stack images. The centers +of the apertures for each object across the 5 filters are fixed so +that the pixels represent the equivalent portions of the same galaxy +for all 5 filters. In this analysis, the best model for each object +is subtracted from the image pixels for all objects excluding the +object in consideration. The 'best model' is \note{TBD}. + +In addition to the raw radial apertures, the stack images are each +convolved with a circular Gaussian with $\sigma$ chosen to yield an +image with a typical FWHM of 6\arcsec. The full set of radial +apertures are again measured on these convolved images. Again, the +best object models are subtracted from the image for objects not being +measured. This subtraction includes the convolution to smooth the +model to the effective FWHM of the convolved image. The entire +procedure is then repeated with a target FWHM of 8\arcsec. + +\note{is the first convolution done with the Alard-Lupton technique?} + \subsection{Forced Photometry : PSFs} \subsection{Forced Photometry : galaxies} -\subsection{Types of Object / PSF models currently available} - -\note{the discussion of the model types needs to be extended} - -\begin{itemize} -\item GAUSS : Pure elliptical Gaussian -\item PGAUSS : polynomial approximation to a Gaussian (PGAUSS) -\item QGAUSS : power law with variable exponential term -\item SGAUSS : -\end{itemize} - -\note{discuss the stability issues with the galaxy model(s)} - \subsection{Output Options} @@ -1395,10 +1518,8 @@ \subsection{Trailed Sources} -\subsection{Fixed / Known-position Sources} - \subsection{Difference Images} The variance map for a difference image must be generated from the two -images use to construct the difference. Otherwise, the low sky level +images used to construct the difference. Otherwise, the low sky level will automatically result in inconsistent interpretation of the variance. @@ -1454,3 +1575,22 @@ \section{Sample Tests} +\begin{verbatim} +Configuration variables affecting the peak detection process: +PEAKS_SMOOTH_SIGMA [2.5] : Gaussian sigma of smoothing kernel, in pixels. +PEAKS_SMOOTH_NSIGMA [2.0] : Gaussian smoothing kernel window size in sigmas. +PEAKS_NSIGMA_LIMIT [20.0] : Detection limit on first pass (sigmas). +PEAKS_NSIGMA_LIMIT_2 [5.0] : Detection limit on faint detection pass (sigmas). +\end{verbatim} + +Recipe parameters which affect the PSF-convolved galaxy model fitting +process: +\begin{verbatim} +EXT_FIT_NSIGMA_CONV [9] : number of sigma +EXT_FIT_ITER +EXT_FIT_MIN_TOL +EXT_FIT_MAX_TOL +LMM_FIT_CHISQ_CONVERGENCE +LMM_FIT_GAIN_FACTOR_MODE +\end{verbatim} + \end{document}