Index: trunk/doc/release.2015/ps1.analysis/analysis.tex =================================================================== --- trunk/doc/release.2015/ps1.analysis/analysis.tex (revision 39945) +++ trunk/doc/release.2015/ps1.analysis/analysis.tex (revision 39946) @@ -85,13 +85,13 @@ \begin{abstract} -Over 3 billion astronomical objects have been detected in the more +Over 3 billion astronomical sources have been detected in the more than 22 million orthogonal transfer CCD images obtained as part of the Pan-STARRS\,1 $3\pi$ survey. Over 85 billion instances of those -objects have been automatically detected and characterized by the +sources have been automatically detected and characterized by the Pan-STARRS Image Processing Pipeline photometry software, \code{psphot}. This fast, automatic, and reliable software was developed for the Pan-STARRS project, but is easily adaptable to images from other telescopes. We describe the analysis of the -astronomical objects by \code{psphot} in general as well as for the +astronomical sources by \code{psphot} in general as well as for the specific case of the 3rd processing version used for the first public release of the Pan-STARRS $3\pi$ survey data. @@ -180,11 +180,11 @@ 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 +measuring the astronomical sources in the image data stream in a timely fashion are a significant challenge. -In order to achieve these ambitious goals, the object detection, +In order to achieve these ambitious goals, the source detection, classification, and measurement process must be both precise and efficient. Not only is it necessary to make a careful measurement of -the flux of individual objects, it is also critical to characterize +the flux of individual sources, it is also critical to characterize the image point-spread-function, and its variations across the field and from image to image. Since comparisons between images must be @@ -192,5 +192,5 @@ astrometry. -A variety of astronomical software packages perform the basic object +A variety of astronomical software packages perform the basic source detection, measurement, and classification tasks needed by the Pan-STARRS IPP. Each of these programs have their own advantages and @@ -211,10 +211,10 @@ automated fashion, does it handle 2D variations well? \citep{1987PASP...99..191S}. -\item Sextractor : pure aperture measurement with rudimentary object +\item Sextractor : pure aperture measurement with rudimentary source subtraction. pro: fast, widely used, easy to automate. con: poor - object separation in crowded regions, PSF-modeling was only in beta, + source separation in crowded regions, PSF-modeling was only in beta, not widely used at the time \citep{sextractor}. -\item galfit : detailed galaxy modeling. not a multi-object PSF +\item galfit : detailed galaxy modeling. not a multi-source PSF analysis tool. con: does not provide a PSF model, not easily automated. very detailed results in very slow processing. only a @@ -233,5 +233,5 @@ re-integrated into the Pan-STARRS pipeline. A new photometry analysis package was developed using lessons learned from the existing -photometry systems. In the process, the object analysis software was +photometry systems. In the process, the source 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 @@ -285,5 +285,5 @@ \begin{itemize} \item {\bf 10 millimagnitude photometric accuracy}. For PSPhot, this - implies that the measured photometry of stellar objects must be + implies that the measured photometry of stellar sources must be substantially better than this 10 mmag since the photometry error per image is combined with an error in the flat-field calibration @@ -299,5 +299,5 @@ astrometric calibration depends on the consistency of the individual measurements. The measurements from PSPhot must be sufficiently - representative of the true object position to enable astrometric + representative of the true source position to enable astrometric calibration at the 10mas level. The error in the individual measurements will be folded together with the errors introduced by @@ -329,6 +329,6 @@ \item {\bf Flexible non-PSF models} PSPhot must be able to represent - PSF-like objects as well as non-PSF sources (e.g., galaxies). It - must be easy to add new object models as interesting representations + PSF-like sources as well as non-PSF sources (e.g., galaxies). It + must be easy to add new source models as interesting representations of sources are invented. @@ -357,18 +357,18 @@ \begin{enumerate} -\item {\bf Image preparation} Load data, characterize the image +\item {\bf Image Preparation} Load data, characterize the image background, load or construct variance and mask images. -\item {\bf Initial object detection} Smooth, find peaks, measure basic +\item {\bf Initial Source Detection} Smooth, find peaks, measure basic properties. -\item {\bf PSF determination} Select PSF candidates, perform model +\item {\bf PSF Determination} Select PSF candidates, perform model fits, build PSF model from fits, select best PSF model class. -\item {\bf Bright object analysis} Fit objects with PSFs, determine - PSF validity, subtract PSF-like objects, fit non-PSF model(s), +\item {\bf Bright Source Analysis} Fit sources with PSFs, determine + PSF validity, subtract PSF-like sources, fit non-PSF model(s), select best model class, subtract model. -\item {\bf Low S/N sources} Detect low-level sources, measure +\item {\bf Faint Source Analysis} Detect low-level sources, measure properties (aperture or PSF) @@ -379,5 +379,5 @@ aperture variations, and background-error corrections. -\item {\bf Output} Write out objects in selected format, write out +\item {\bf Output} Write out sources in selected format, write out difference image, variance image, etc, as selected. \end{enumerate} @@ -389,8 +389,10 @@ case the PSF modeling stage can be skipped. +{\bf A note on nomenclature:} + \subsection{Image Preparation} The first step is to prepare the image for detection of the -astronomical objects. We need three separate images: the measured +astronomical sources. We need three separate images: the measured flux (signal image), the corresponding variance image, and a mask defining which pixels are valid and which should be ignored. The @@ -405,6 +407,6 @@ be constructed automatically by PSPhot. -The mask is represented as 16-bit integer image in which a value of 0 -represents a valid pixel. Each of the 16 bits define different +The mask is represented as a 16-bit integer image in which a value of +0 represents a valid pixel. Each of the 16 bits define different reasons a pixel should be ignored. This allows us to optionally respect or ignore the mask depending on the circumstance. For @@ -413,19 +415,26 @@ saturated pixel. In addition, the mask pixels are used to define the pixels available during a model fit, and which should be ignored for -that specific fit. The initial mask, if not supplied by the user, is -constructed by default from the image by applying three rules: 1) -Pixels which are above a specified saturation level are marked as -saturated (configuration keyword: \code{SATURATE}). 2) Pixels which -are below a user-defined value are considered unresponsive and masked -as dead. 3) Pixels which lie outside of a user-defined coordinate -window are considered non-data pixels (eg, overscan) and are marked as -invalid. The valid window is defined by the configuration variables -\code{XMIN}, \code{XMAX}, \code{YMIN}, \code{YMAX}. - -PSPhot (and other IPP) functions understand two types of masked -pixels: ``bad'' and ``suspect''. Bad pixels are those which should -not be used in any operations, while suspect pixels are those for -which the reported signal may be contaminated or biased, but may be -useable in some contexts. For example, a pixel with poor charge +that specific fit (\code{MARK = 0x8000}). The initial mask, if not +supplied by the user, is constructed by default from the image by +applying three rules: 1) Pixels which are above a specified saturation +level are marked as saturated. The level is specified by the camera +format keyword \code{CELL.SATURATION}, which may specify a value or +define a header keyword which in turn specifies the value in the image +header. In the case of PS1 PV3, the header keyword \code{MAXLIN} +specifies the saturation level for each chip. \note{refer to detrend + paper here? what are GPC1 saturation levels?}. 2) Pixels which are +below a user-defined value are considered unresponsive and masked as +dead. (camera format keyword \code{CELL.BAD} = 0 for PS1 PV3). 3) +Pixels which lie outside of a user-defined coordinate window are +considered non-data pixels (eg, overscan) and are marked as invalid. +(psphot recipe keywords \code{XMIN}, \code{XMAX}, \code{YMIN}, +\code{YMAX}, all set to 0 for PS1 PV3 -- invalid pixels were specified +for PS1 PV3 with a supplied mask image, see \cite{waters2017}. + +The library functions used by \code{psphot} understand two types of +masked pixels: ``bad'' and ``suspect''. Bad pixels are those which +should not be used in any operations, while suspect pixels are those +for which the reported signal may be contaminated or biased, but may +be useable in some contexts. For example, a pixel with poor charge transfer efficiency is likely to be too untrustworthy to use in any circumstance, while a pixel in which persistence ghosts have been @@ -465,16 +474,18 @@ \end{table*} -The variance image, if not supplied is constructed by default from the -flux image using the configuration supplied values of \code{GAIN} and -\code{READ_NOISE} to calculate the appropriate Poisson statistics for -each pixel. In this case, the image is assumed to represent the -readout from a single detector, with well-defined gain and read noise -characteristics. This assumption is not always valid. For example, -if the input flux image is the result of an image stack with a -variable number of input measurements per pixel (due to masking and -dithering), the variance cannot be calculated from the signal image -alone. It is necessary in such a case to supply a variance image which -accurately represents the variance as a function of position in the -image. +The variance image, if not supplied, is constructed by default from +the flux image using the configuration supplied gain and read noise +values to calculate the appropriate Poisson statistics for each pixel. +The parameters are determined based on the camera format keywords +\code{CELL.GAIN} and \code{CELL.READNOISE}, which in the case of PS1 +PV3 refer to the header keywords \code{GAIN} and \code{RDNOISE}. In +this case, the image is assumed to represent the readout from a single +detector, with well-defined gain and read noise characteristics. This +assumption is not always valid. For example, if the input flux image +is the result of an image stack with a variable number of input +measurements per pixel (due to masking and dithering), the variance +cannot be calculated from the signal image alone. It is necessary in +such a case to supply a variance image which accurately represents the +variance as a function of position in the image. Some image processing steps introduce cross-correlation between pixel @@ -489,25 +500,61 @@ covariance image is prohibitive. +\note{describe the way we handle covariance} + 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. 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. See also the discussion in -\cite{waters2017}. - -\subsection{Initial Object Detection} +with a spacing defined by the psphot recipe values +\code{BACKGROUND.XBIN, BACKGROUND.YBIN}, set to 400 pixels for PS1 +PV3. Superpixels of size \code{BACKGROUND.XSAMPLE, + BACKGROUND.YSAMPLE} ($2 \times 2$ for PS1 PV3) times larger than +this spacing are used to measure the local background for each +background grid point, thus over-sampling the background spatial +variations. In the interest of speed, a subset of \code{IMSTATS_NPIX} +(10,000 for PS1 PV3) randomly selected {\em unmasked} pixels in these +regions are used to determine the background. The background value +for each superpixel is determined by fitting a Gaussian distribution +to the histogram of pixels values. + +If the image were empty of stars and only contained flux from a +uniform background sky, we would expect the distribution to be Poisson +distributed, and in general in a high-enough signal range to be +essentially Gaussian. We fit a symmetric Gaussian to all histogram +bins within 15\% of the peak bin value to determine the mean and +standard deviation values for the background. + +If, however, the sky is not empty of stars or other sources, and we +have correctly masked the large majority of non-responsive pixels, +then we expect the flux distribution of the pixels to be asymmetric +with a Gaussian core representing the sky and a tail to the high end +representing the pixels with astronomical source flux contributions. +We would like to determine the mean of the underlying Gaussian without +suffering bias from the stellar flux. We thus perform a second +Gaussian fit using an asymmetric subset of the histogram pixels, +fitting those histogram bins which are left of the peak but above 25\% of +the peak value, or right of the peak but above 50\% of the peak +value. + +If the fit to the asymmetric lower fraction of the curve is less than +the symmetric fit, but greater than the above lower-bound of the full +symmetric fit, then the lower fraction value is kept as the true mean +sky value for this superpixel. + +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 source in the +output catalog. See also the discussion in \cite{waters2017}. + +\note{give examples with simulations and show examples of over-subtraction} + +\subsection{Initial Source Detection} \subsubsection{Peak Detection} \label{sec:peaks} -The objects are initially detected by finding the location of local +\note{add a ref to the Kaiser paper} + +The sources 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. The @@ -516,8 +563,8 @@ 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. +keyword: \code{PEAKS_NSIGMA_LIMIT} = 20.0 for PS1 PV3). A maximum of +\code{PEAKS_NMAX} (5000 of PS1 PV3) 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 @@ -529,5 +576,5 @@ any of the other 8 pixels is kept if the pixel $X$ and $Y$ coordinates are greater than or equal to the other equal value pixels. This -simple rule set means that a flat-topped region will maintain peaks at +simple rule set means that a flat-topped region will result peaks at the maximum $X$ and $Y$ corners of the region. @@ -585,4 +632,9 @@ \end{eqnarray} +The resulting peak position, ($x_{min}, y_{min}$), is used as the +default starting coordinate for the source. Later in the +\code{psphot} analysis, improved measurements of the source positions +are calculated as discussed below. + \begin{figure}[htbp] \begin{center} @@ -601,10 +653,11 @@ 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). +the smoothed significance image above the detection threshold +(\code{PEAKS_NSIGMA_LIMIT}). 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 \code{FOOTPRINT_GROW_RADIUS} (= 3 +pixels for PS1 PV3). Peaks are assigned to the footprints in which +they are contained (note by construction all peaks must be located in +a footprint since the peaks must be above the detection threshold). For any peak which is not the brightest peak in that footprint it is @@ -613,16 +666,17 @@ {\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 (see +\code{FOOTPRINT_CULL_NSIGMA_DELTA} (4.0 for PS1 PV3) sigmas below the +peak of interest, the peak is considered to be {\em locally + insignificant} and removed from the list of possible detections (see Figure~\ref{fig:peaks}). In the vicinity of a saturated star, the -rule is somewhat more agressive as the flat-topped or structured +rule is somewhat more aggressive 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} +the proximity to saturation. Sources for which the peak is greater +than 50\% of the saturation value 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} \label{sec:moments} @@ -645,32 +699,33 @@ 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. +properties of the sources related to the first, second, and higher +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, -$\sigma_{\rm PSF}$. This -window function reduces the noise of the measurement of the first and -second moments by suppressing the noisy pixels at high radial distance -as well as by reducing the contaminating effects of neighboring stars. -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_{\rm PSF}$ which is modified by the effect of -the window function. In addition, depending on the size of the window -function compared to the true PSF size, the measured value of the PSF -size, $\sigma_{\rm PSF}$, will be biased high or low depending on the -signal-to-noise of the object. +a ``window function'', down-weighting the pixels by a Gaussian, +centered on the object, with size $\sigma_W$ chosen to be large +compared to the PSF size, $\sigma_{\rm PSF}$. This window function +reduces the noise of the measurement of the moments by suppressing the +noisy pixels at high radial distance as well as by reducing the +contaminating effects of neighboring stars. The choice of $\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 the PSF size, +$\sigma^{\prime}_{\rm PSF}$ which different from the true value due to +the effect of the window function. The measured value of the PSF size +will be biased high or low depending on both the signal-to-noise of +the source and the size of the window function compared to the true +PSF size. These effects are illustrated in Figure~\ref{fig:moments.window} using simulated data. An image was generated with a PSF model matching the -radial profile of the PS1 PSF model with a FWHM of 1.4 arcseconds. As -the window function $\sigma_W$ is increased, the measured FWHM for the -bright simulated stars rises to meet the truth value. For small -values of $\sigma_W$, fainter stars are biased to low measured values -of the FWHM. For large values of $\sigma_W$, the faint stars are -biased to higher values and the scatter increases. +radial profile of the PS1 PSF model with $\sigma_{\rm PSF}$ +corresponding to a FWHM of 1.4 arcseconds. As the window function +$\sigma_W$ is increased, the measured FWHM for the bright simulated +stars rises to meet the truth value. For small values of $\sigma_W$, +fainter stars are biased to low measured values of the FWHM. For +large values of $\sigma_W$, the faint stars are biased to higher +values and the scatter increases. In a real image, we do not know the true value of the PSF size. If we @@ -681,41 +736,46 @@ artifacts) and (2) the brighter stars are themselves subject to additional biases due to saturation and other non-linear effects -(c.f., ``the Brighter-Fatter'' effect, REF). To make a robust -choice for the window function $sigma_w$, we choose a value -such that the measured value of $\sigma_{\rm PSF}$ is 65\% of -$\sigma_w$. The resulting second moment values are biased somewhat -low (\approx 75\% of the truth value for the PS1 PSF profile), but are -relatively unbiased as a function of brightness. - -To choose the value of $\sigma_W$, we try values of (1, 2, 3, 4.5, 6, -9, 12, 18) pixels $\approx$ (0.26, 0.51, 0.77, 1.15, 1.54, 2.3, 3.1, -4.6) arcseconds. 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$, we 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. -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. +(c.f., ``the Brighter-Fatter'' effect, \note{REF}). To make a robust +choice for $\sigma_w$, we choose a value such that the measured value +of $\sigma^{\prime}_{\rm PSF}$ is 65\% of $\sigma_w$. The resulting second +moment values are biased somewhat low (\approx 75\% of the truth value +for the PS1 PSF profile), but are relatively unbiased as a function of +brightness. + +To choose the value of $\sigma_W$, we try a sequence of values +spanning a range guaranateed to contain any reasonable seeing values. +The values are specified in the \code{psphot} recipe as +\code{PSF.SIGMA.VALUES} and have the following values for PS1 PV3: (1, +2, 3, 4.5, 6, 9, 12, 18) pixels $\approx$ (0.26, 0.51, 0.77, 1.15, +1.54, 2.3, 3.1, 4.6) arcseconds. For each of these $\sigma_W$ values, +we then select candidate PSF stars based on the distribution of the +measured $\sigma^{\prime}_{\rm PSF}$ in the two principal directions: +$\sigma_{x,x}$ and $\sigma_{y,y}$ (see +Section~\ref{sec:psf.source.selection}, below). For each test value +of $\sigma_w$, we determine the ratio $\rho_\sigma = +\frac{\sigma_{x} + \sigma{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 $\rho_\sigma$ is expected to be 0.65. We use +an aperture with a radius of 4$\sigma_w$ to select the pixels for the +measurement of the moments. + +Once $\sigma_w$ 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 +x_0 & = & \frac{1}{S} \sum_i w_i (f_i - s_i)x_i \\ +y_0 & = & \frac{1}{S} \sum_i w_i (f_i - s_i)y_i \\ +M_{xx} & = & \frac{1}{S} \sum_i w_i (f_i - s_i)(x_i - x_0)^2 \\ +M_{xy} & = & \frac{1}{S} \sum_i w_i (f_i - s_i)(x_i - x_0)(y_i - y_0) \\ +M_{yy} & = & \frac{1}{S} \sum_i w_i (f_i - s_i)(y_i - y_0)^2 \\ +M_{xxx} & = & \frac{1}{S} \sum_i \frac{w_i}{r_i} (f_i - s_i)(x_i - x_0)^3 \\ +M_{xxy} & = & \frac{1}{S} \sum_i \frac{w_i}{r_i} (f_i - s_i)(x_i - x_0)^2(y_i - y_0) \\ +M_{xyy} & = & \frac{1}{S} \sum_i \frac{w_i}{r_i} (f_i - s_i)(x_i - x_0)(y_i - y_0)^2 \\ +M_{yyy} & = & \frac{1}{S} \sum_i \frac{w_i}{r_i} (f_i - s_i)(y_i - y_0)^3 \\ +M_{xxxx} & = & \frac{1}{S} \sum_i \frac{w_i}{r^2_i} (f_i - s_i)(x_i - x_0)^4 \\ +M_{xxxy} & = & \frac{1}{S} \sum_i \frac{w_i}{r^2_i} (f_i - s_i)(x_i - x_0)^3(y_i - y_0) \\ +M_{xxyy} & = & \frac{1}{S} \sum_i \frac{w_i}{r^2_i} (f_i - s_i)(x_i - x_0)^2(y_i - y_0)^2 \\ +M_{xyyy} & = & \frac{1}{S} \sum_i \frac{w_i}{r^2_i} (f_i - s_i)(y_i - y_0)(y_i - y_0)^3 \\ +M_{yyyy} & = & \frac{1}{S} \sum_i \frac{w_i}{r^2_i} (f_i - s_i)(y_i - y_0)^4 \end{eqnarray} where $f_i$ is the flux in a pixel; $s_i$ is the local sky value for @@ -723,19 +783,19 @@ $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, +pixel, $\sqrt{(x_i - x_0)^2 + (y_i - y_0)^2}$; The sums are performed +over all (unmasked) 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: +If the measured centroid coordinates ($x_0, y_0$) differ from the peak +coordinates be a large amount (1.5$\sigma_w$), then the peak is +identified as being of poor quality (\code{infoFlag} bit +\code{MOMENTS_FAILURE}) and is skipped in further analyses. In such +as case, 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, the 1st and half-radial moments, +$M_r$ and $M_h$ as defined below, are calculated: \begin{eqnarray} M_r & = & \frac{1}{S} \sum_i (f_i - s_i)r_i \\ @@ -745,28 +805,39 @@ these moments. -The Kron radius \citep{1980ApJS...43..305K} 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. +With the first radial moment, we can calculate a preliminary Kron +radius and magnitude. The Kron radius \citep{1980ApJS...43..305K} 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\sigma_w$ if that +cannot be determined. $R_{\rm max}$ is set to the size of the moments +aperture, $4\sigma_w$. At this stage, the measurement of the Kron +parameters are preliminary since the aperture has been chosen as a +fixed size relative to the size of the PSF. At a later stage, +higher-quality Kron parameters appropriate to galaxies are measured +with more care paid to the exact aperture used +(Section~\ref{sec:kron.mags}). + +% $\sigma_w$ is saved as MOMENTS_GAUSS_SIGMA +% the aperture radius is saved as PSF_MOMENTS_RADIUS \subsection{PSF Determination} -\subsubsection{PSF Model vs Object Model} - -PSPhot uses an analytical model to represent the shape and flux of an -object. An important concept within the PSPhot code is the -distinction between a model which describes an object on an image and -a model with describes the point-spread-function (PSF) across an -image. - -Any object in an image may be represented by some analytical model, +\subsubsection{PSF Model vs Source Model} + +The PSF model used by \code{psphot} consists of an analytical function +combined with a pixelized representation of the residual differences +between the analytical model and the true PSF. Both the shape +parameters of the analytical model and the pixelized residual +differences are allowed to vary in two dimensions across the images. + +Within \code{psphot}, several analytical models may be used to +describe the PSF, but all share a few common characteristics. + +Any source in an image may be represented by some analytical model, for example, a 2-D elliptical Gaussian: \begin{eqnarray} @@ -776,24 +847,24 @@ y & = & y_{\rm ccd} - y_o \end{eqnarray} -The object model will have a variety of model parameters, in this case +The source model will have a variety of model parameters, in this case the centroid coordinates ($x_o, y_o$), the elliptical shape parameters ($\sigma_x, \sigma_y, \sigma_{\rm xy}$), the model normalization ($I_o$) and the local value of the background ($S$). A specific -object will have a particular set of values for these different +source will have a particular set of values for these different parameters. 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 +unresolved sources in the image. In a typical image, the shape of 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 +the position of the source on the image. Note that the source 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 +source to source. 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$). -Thus these parameters are each a function of the object centroid +Thus these parameters are each a function of the source centroid coordinates: \begin{eqnarray} @@ -806,5 +877,5 @@ configuration. The first option is to use a 2-D polynomial which is fitted to the measured parameter values across the image. The second -option is to use a grid of values which are measured for objects +option is to use a grid of values which are measured for sources within a subregion of the image. In the latter case, the value at a specific coordinate in the image is determined by interpolation @@ -822,7 +893,7 @@ % XXX discuss the improvements in the astrometric modeling PV1 - PV3 -PSPhot uses a single structure to represent the object model and -another structure to represent the PSF model. The object model -structure consists of the collection of measured object model +PSPhot uses a single structure to represent the source model and +another structure to represent the PSF model. The source model +structure consists of the collection of measured source model parameters, carried as a \code{psLib} vector (\code{psVector}) along with an equal-length vector with the parameter errors. The structure @@ -834,21 +905,21 @@ The PSPhot representation of the PSF consists of an array of -polynomials, each representing the variation in the object model PSF +polynomials, each representing the variation in the source model PSF parameters (\code{psArray} of \code{psPolynomial2D}). The PSF model structure also includes the same integer used to identify which model corresponds to particular instance of the PSF. At the moment, the number of PSF parameters is a fixed number (4) fewer than the number -of parameters of the corresponding object model. For example, the -elliptical Gaussian model uses 7 parameters to represent the object and +of parameters of the corresponding source model. For example, the +elliptical Gaussian model uses 7 parameters to represent the source and 3 for the PSF model. -PSPhot is written so that the object detection, measurement, and +PSPhot is written so that the source detection, measurement, and classification code does not depend on the specific form of the -available object model functions. Access to the characteristics of +available source model functions. Access to the characteristics of the models is provided through a simple function abstraction method. Throughout PSPhot, there are many places where it is necessary for the -code to refer to an aspect of the object or PSF model. Often, these +code to refer to an aspect of the source or PSF model. Often, these quantities are needed deep within other parts of the code. For -example, when attempting to fit the pixel flux values for an object, +example, when attempting to fit the pixel flux values for a source, it is necessary to generate a guess for the model parameters. Or, in order to limit the domain of the fit, it is necessary to determine an @@ -872,24 +943,25 @@ When a new model is provided to PSPhot, it is not necessary to specify -the intended use of the object model function (ie, PSF-like object, +the intended use of the source model function (ie, PSF-like source, galaxy, comet, etc). Any model can be used for the PSF model, or to -describe the flux distributions of the non-PSF objects. The code -currently uses a fixed translation between the object model parameters +describe the flux distributions of the non-PSF sources. The code +currently uses a fixed translation between the source model parameters and the PSF model parameters. It also defines a specific order for the 4 independent parameters. -\subsubsection{PSF Candidate Object Selection} +\subsubsection{Candidate PSF Source Selection} +\label{sec:psf.source.selection} The first stage of determining the PSF model for an image is to -identify a collection of objects in the image which are {\em likely} -to be PSF-like. PSPhot uses the object moments to make the initial -guess at a collection of PSF-like objects. At this point, the program -has measured the second order moments for all objects identified by +identify a collection of sources in the image which are {\em likely} +to be PSF-like. PSPhot uses the source moments to make the initial +guess at a collection of PSF-like sources. At this point, the program +has measured the second order moments for all sources identified by their peaks, as well as an approximate signal-to-noise ratio. All -objects with a S/N ratio greater than a user-defined parameter +sources with a S/N ratio greater than a user-defined parameter (\code{PSF_SHAPE_NSIGMA} = 20.0) are selected by PSPhot, though -objects which have more than a certain number of saturated pixels are +sources which have more than a certain number of saturated pixels are excluded at this stage. PSPhot then examines the 2-D plane of -$\sigma_x, \sigma_y$ in search of a concentrated clump of objects (see +$\sigma_x, \sigma_y$ in search of a concentrated clump of sources (see Figure~\ref{fig:moment.class}). To do this, it constructs an artificial image with pixels representing @@ -901,15 +973,15 @@ is then examined to find a peak which has a significance greater than XXX. Unless the image is extremely sparse, such a peak will be -well-defined and should represent the objects which are all very -similar in shape. Other objects in the image will tend to land in +well-defined and should represent the sources which are all very +similar in shape. Other sources in the image will tend to land in very different locations, failing to produce a single peak. To avoid -detecting a peak from the unresolved cosmic rays, objects which have +detecting a peak from the unresolved cosmic rays, sources which have second-moments very close to 0 are ignored. The only danger is if the -PSF is very small and too many of these objects are rejected as cosmic +PSF is very small and too many of these sources are rejected as cosmic rays. Once a peak has been detected in this plane, the centroid and second -moments of this peak are measured. All objects which land within XXX -$\sigma$ of this centroid are selected as likely PSF-like objects in +moments of this peak are measured. All sources which land within XXX +$\sigma$ of this centroid are selected as likely PSF-like sources in the image. @@ -927,9 +999,9 @@ \end{figure} -\subsubsection{PSF Candidate Object Model Fits} +\subsubsection{Candidate PSF Source Model Fits} % \note{link to psLibADD} -All candidate PSF objects are then fitted with the selected object +All candidate PSF sources are then fitted with the selected source model, allowing all of the parameters (PSF and independent) to vary in the fit. PSPhot uses the Levenberg-Marquardt minimization technique @@ -938,8 +1010,8 @@ starting 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. Any objects which fail to +shape parameters for the PSF models. Any sources which fail to converge in the fit are flagged as invalid. -For the resulting collection of object model parameters, the +For the resulting collection of source 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 a @@ -948,14 +1020,14 @@ 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 +without being excessively biased by individual sources which fail +drastically. Sources 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 -model to specify the dependent model parameter values for each object. +All of the PSF-candidate sources are then re-fitted using the PSF +model to specify the dependent model parameter values for each source. For example, in the case of the elliptical Gaussian model, the shape -parameters ($\sigma_x, \sigma_y, \sigma_{xy}$) for each object are -set by the coordinates of the object centroid and fixed (not allowed +parameters ($\sigma_x, \sigma_y, \sigma_{xy}$) for each source are +set by the coordinates of the source centroid and fixed (not allowed to vary) in the fitting procedure. The resulting fitted models are then used to determine a metric which tests the quality of the PSF @@ -964,9 +1036,9 @@ The metric used by PSPhot to assess the PSF model is the scatter in the differences between the aperture and fit magnitudes for the PSF -objects. The difference between the aperture and fit magnitudes ({\em +sources. The difference between the aperture and fit magnitudes ({\em ApResid}) is a critical parameter for any PSF modeling software which uses an analytical model to represent the flux distribution of the -objects in an image. An approximate correction is measured here, with -a more detailed correction measured after all object analysis is +sources in an image. An approximate correction is measured here, with +a more detailed correction measured after all source analysis is performed. The PSF model with the best consistency of the aperture correction is judged to be the best model. @@ -974,14 +1046,44 @@ \subsection{Bright Source Analysis} -%% \subsubsection{Very Bright Stars} -%% -%% The PSF modeling code fails to fit the wings of highly saturated stars -%% if the core of the star is too contaminated by saturated pixels. For -%% stars with estimated instrumental magnitudes brighter than XXX, we fit -%% and subtract a radial profile modeled with a spline (?). +Once a PSF model has been determined, the brighter sources in the +image may be analysed in detail. The goals in this stage are (1) to +determine the fluxes and positions of the bright stellar sources with +high precision appropriate to their high signal-to-noise and (2) to +characterize the bright source flux profiles sufficiently well that +they may be subtracted from the image to allow for the clean detection +of the fainter sources. Note that as the analysis proceeds, there are +several stages in which the 2D flux models for all sources are +subtracted from the image, and individual sources are replaced in the +image for a particular analysis step and then removed again. + +In order to allow for multiple threads to process a single image, the +pixels in an image are divided into a grid of superpixels (see +Figure~\ref{fig:threadgrid}). The superpixels are assigned to one of +four groups, as illustrated, so that each superpixel in a group is +well separated from the other superpixels of that group. The analysis +of the image proceeds in 4 steps, one for each of these groups. Each +of the superpixels in the first group is assigned to a single thread +until all threads are assigned. A single thread is responsible for +the analysis of sources which land within their current superpixel, as +determined by the centroid coordinates. As the threads complete their +analysis, they are assigned the next unfinished superpixel in the +active group. When all superpixels in one group have been processed, +then the superpixels in the next group can start. This strategy +allows the threading to process sources which may be extended without +the danger that two threads are actively touching the same pixels. +For the PV3 analysis, 4 threads were used for most processing tasks. + +\subsubsection{Very Bright Stars} + +The PSF modeling code fails to fit the wings of highly saturated stars +if the core of the star is too contaminated by saturated pixels. For +stars with estimated instrumental magnitudes brighter than XXX, we fit +and subtract a radial profile modeled with a spline (?). + +\note{more here} \subsubsection{Fast Ensemble PSF Fitting} -Before the detailed analysis of the objects is performed, it is +Before the detailed analysis of the sources is performed, it is 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 @@ -1019,53 +1121,53 @@ achieve a good convergence. -Once a solution set for $A_i$ is found, all of the objects are +Once a solution set for $A_i$ is found, all of the sources are subtracted from the by applying these values to the unit-flux PSF. -\subsubsection{PSF Model applied to detected objects} +\subsubsection{Full PSF Model Fitting} % \note{review the discussion below} Once a PSF model has been selected for an image, PSPhot attempts to -fit all of the detected objects, above a user-defined signal-to-noise +fit all of the detected sources, above a user-defined signal-to-noise ratio with the PSF model. For these fits, the dependent parameters -are fixed by the PSF model and only the 4 independent object model +are fixed by the PSF model and only the 4 independent source model parameters are allowed to vary in the fit. PSPhot again uses Levenberg-Marquardt minimization for the non-linear fitting. The -objects are fitted in their S/N order, starting with the brightest and +sources are fitted in their S/N order, starting with the brightest and working down to the user-specified limit. -Once a solution has been achieved for an object, PSPhot attempts to -judge the quality of the PSF model as a representation of the object +Once a solution has been achieved for a source, PSPhot attempts to +judge the quality of the PSF model as a representation of the source shape. To do this, it calculates the next step of the minimization {\em allowing 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 -large if the PSF is not a good representation of the object flux. The +should be very small if the source is well represented by the PSF, but +large if the PSF is not a good representation of the source flux. The model quality is judged by the change in the two shape parameters -which represent the 2D size of the object. For the case of the +which represent the 2D size of the source. For the case of the elliptical Gaussian, these two parameters are $\sigma_x$ and $\sigma_y$. For a generic model, the shape parameters may be defined differently, but there should always be two parameters which scale the -object size in two dimensions. Currently, PSPhot requires the two +source 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 +will be a function of the signal-to-noise of the source in question and the value of the shape parameter itself. The expected standard 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 +the source is well-represented by the PSF, then the shape parameter values should be close to their minimization value. We can thus ask, -for each object, given the measured amplitude of the Gauss-Newton +for each source, given the measured amplitude of the Gauss-Newton step, how many standard deviations from the expected value (of 0.0) is -this particular value? Objects for which the variation in the shape +this particular value? Sources for which the variation in the shape parameters is a large positive number of standard deviations are likely to be better represented by a larger flux distribution than the -PSF (eg, a Galaxy or Comet, etc). Objects for which the variation in +PSF (eg, a Galaxy or Comet, etc). Sources for which the variation in the shape parameters is a large negative number of standard deviations are likely to be better represented by a smaller flux distribution than the PSF (ie, a cosmic ray or other defect). A user-defined number of standard deviations is used to select these two cases, and -to flag the object as a likely galaxy (really meaning 'extended') or +to flag the source as a likely galaxy (really meaning 'extended') or as a likely defect. @@ -1082,9 +1184,9 @@ converge on a fit with very low or negative peak flux / flux normalization. PSPhot will flag any non-convergent PSF fit and any -object with PSF S/N ratio lower than a user-defined cutoff. It is +source with PSF S/N ratio lower than a user-defined cutoff. It is also useful to identify very poor fits by setting a maximum Chi-Square -cutoff for objects. - -As the objects are fitted to the PSF model, those which survive the +cutoff for sources. + +As the sources are fitted to the PSF model, those which survive the exclusion stage are subtracted from the image. The subtraction process modifies the image pixels (removing the fitted flux, though @@ -1100,5 +1202,5 @@ Sources which are blended with other sources are fitted together as a set of -PSFs. A single multi-object fit is performed on all blended peaks. +PSFs. A single multi-source fit is performed on all blended peaks. The resulting fits are evaluated independently and any which are determined to be PSFs are subtracted from the image. @@ -1107,5 +1209,5 @@ Sources which are judged to be non-PSF-like are confronted with two -possible alternative choices. First, the object is fitted with a +possible alternative choices. First, the source is fitted with a double-source model. In this pass, the assumption is made that there are two neighboring sources, but the peaks are blended together, or @@ -1125,8 +1227,8 @@ 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 +of the source. If the source 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 +source model in special cases). If the source is small compared to a PSF, it is considered to be a {\em cosmic ray} and masked. @@ -1134,13 +1236,13 @@ 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 +and Kron magnitudes, is calculated for each source. 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. +the source 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 width of the source in the minor axis direction. The second-moment in the minor axis direction is calculated from $M_{xx}, M_{xy}, M_{yy}$ as follows: @@ -1149,11 +1251,11 @@ \] 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 +magnitude is $< -5.5$, then the source is identified as a cosmic ray and the associated pixels are masked. -\subsubsection{Non-PSF Objects} - -Once every object (above the S/N cutoff) has been confronted with the -PSF model, the objects which are thought to be galaxies (extended) can +\subsubsection{Non-PSF Sources} + +Once every source (above the S/N cutoff) has been confronted with the +PSF model, the sources which are thought to be galaxies (extended) can now be fit with appropriate models for the galaxies (or other likely extended shapes). Again, the fitting stage starts with the brightest @@ -1163,15 +1265,15 @@ PSPhot will use the user-selected galaxy model to attempt the galaxy model fits. In the configuration system, the keyword \code{GAL_MODEL} -is set to the model of interest. All suspected extended objects are +is set to the model of interest. All suspected extended sources are fitted with the model, allowing all of the parameters to float. The initial parameter guesses are critical here to achieving convergence on the model fits in a reasonable time. The moments and the pixel flux distribution are used to make the initial parameter guess. Many -of the object parameters can be accurately guessed from the first and +of the source parameters can be accurately guessed from the first and second moments. The power-law slope can be guessed by measuring the isophotal level at two elliptical radii and comparing the ratio to that expected. -For each of the galaxy models (in fact for all object models), a +For each of the galaxy models (in fact for all source models), a function is defined which examines the fit results and determines if the fit can be consider as a success or a failure. The exact criteria @@ -1184,9 +1286,9 @@ All galaxy model fits which are successful are then subtracted from the image as is done for the successful PSF model fits. Of course, -the background flux is retained, with the result that only the object +the background flux is retained, with the result that only the source is subtracted from the image. Again, the variance image is (currently) not modified. -\subsection{Faint Sources} +\subsection{Faint Source Analysis} After a first pass through the image, in which the brighter sources @@ -1194,142 +1296,25 @@ subtracted, PSPhot optionally begins a second pass at the image. In this stage, the new peaks are detected on the image with the bright -objects subtracted. In this pass, the peak detection process uses the +sources subtracted. In this pass, the peak detection process uses the variance image to test the validity of the individual peaks. All peaks with a significance greater than a user-defined minimum threshold are -accepted as objects of potential interest. - -The objects which are measured in this faint-object stage are clearly -low significance detections. The PV3 threshold for the bright object +accepted as sources of potential interest. + +The sources which are measured in this faint-source stage are clearly +low significance detections. The PV3 threshold for the bright source analysis is a signal-to-noise of 20. The lower limit cutoff for the -faint object analysis in PV3 is a signal-to-noise of 5.0. Objects -detected in the faint object stage are fitted with the PSF model using +faint source analysis in PV3 is a signal-to-noise of 5.0. Sources +detected in the faint source stage are fitted with the PSF model using the linear, ensemble fitting process. -\subsection{Aperture Correction Measurement} - -The important concept here is that an analytical model will always -fail to describe the flux of the objects at some level. In the end, -all astronomical photometry is in some sense a relative measurement -between two images. Whether the goal is calibration of a science -image taken at one location to a standard star image at another -location, or the goal is simply the repetitive photometry of the same -star at the same location in the image, it is always necessary to -compare the photometry between two images. If this measurement is to -be consistent, then the measurement must represent the flux of the -stars in the same way regardless of the conditions under which the -images were taken, at least within some range of normal image -conditions. So, for example, two images with different image quality, -or with different tracking and focus errors, will have different PSF -models. Since an analytical model will always fail to represent the -flux of the star at some level, the measured flux of the same object -in the two images will be different (even assuming all other -atmospheric and instrumental effects have been corrected). The -amplitude of the error will by determined by how inconsistently the -models represent the actual object flux. For example, if the first -image PSF model flux is consistently 10\% too low and the second is 5\% -too high, then the comparison between the two images will be in error -by 15\%. - -Aperture photometry avoids these problems, by trading for other -difficulties. In aperture photometry, if a large enough aperture is -chosen, the amount of flux which is lost will be a small fraction of -the total object flux. Even more importantly, as the image conditions -change, the amount lost will change by an even smaller fraction, at -least for a large aperture. This can be seen by the fact that the -dominant variations in the image quality are in the focus, tracking -and seeing. All of these errors initially affect the cores of the -stellar images, rather than the wide wings. The wide wings are -largely dominated by scattering in the optics and scattering in the -atmosphere. The amplitude and distribution of these two scattering -functions do not change significantly or quickly for a single -telescope and site. - -The difficulty for aperture photometry is the need to make an accurate -measurement of the local background for each object. As the aperture -grows, errors in the measurement of the sky flux start to become -dominant. If the aperture is too small, then variation in the image -quality are dominant. The brighter is the object, the smaller is the -error introduced by the large size of the aperture. However, the -number of very bright stars is limited in any image, and of course the -brighter stars are more likely to suffer from non-linearity or -saturation. PSPhot measures the aperture correction ({\em ApResid}) -for every PSF candidate object and applies this correction to the PSF -model photometry. - -% 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 -%%% candidate object, then calculates the trend of this correction as a -%%% function of the magnitude. This trend is fitted with a line. The -%%% resulting function can be used to determine the effective aperture -%%% correction for an infinite flux object and the average bias inherent -%%% in the sky measurement for the image. The scatter of the -%%% PSF-candidate object measurements about this trend is a measure of how -%%% well we can measure photometry from the image by applying the specific -%%% PSF model. The slope of this trend is a measure of the bias in the -%%% local sky measurment for each object. In principal, the measured sky -%%% levels could be modified by this bias. More generally, the measured -%%% bias in a collection of images could be used to improve the model -%%% fitting or sky fitting portion of the software the remove the bias -%%% term. - -PSPhot allows a collection of PSF model functions to be tried on all -PSF candidate objects. For each model test, the above corrected -ApResid scatter is measured. The PSF model function with the smallest -value for the ApResid scatter is then used by PSPhot as the best PSF -model for this image. The number of models to be tested is specified -by the configuration keyword \code{PSF_MODEL_N}. The configuration -variables \code{PSF_MODEL_0}, \code{PSF_MODEL_1}, through -\code{PSF_MODEL_N - 1} specify the names of the models which should be -tested. - -Several likely PSF model classes are available within \code{psphot}: -\begin{itemize} -\item Gaussian : $f = I_0 e^{-z}$ -\item Pseudo-Gaussian : $f = I_0 (1 + z + \frac{1}{2} z^2 + \frac{1}{6} z^3)^{-1}$ \code{[PGAUSS]} -\item Variable Power-Law : $f = I_0 (1 + z + z^{\alpha})^{-1}$ \code{[RGAUSS]} -\item Steep Power-Law : $f = I_0 (1 + \kappa z + z^{2.25})^{-1}$ \code{[QGAUSS]} -\item PS1 Power-Law : $f = I_0 (1 + \kappa z + z^{1.67})^{-1}$ \code{[PS1_V1]} -\end{itemize} -where $z \propto r^2$ ($z = \frac{x^2}{2\sigma_x^2} + -\frac{y^2}{2\sigma_y^2} + \sigma_{\rm xy} x y $). The Pseudo-Gaussian -is a Taylor expansion of the Gaussian and is used by Dophot -\citep{1993PASP..105.1342S}. The latter profiles are similar to the -Moffat profile form \citep{1969AA.....3..455M,1983AA...126..278B}, -with small differences. For the PS1 GPC1 analysis, we used the -\code{PS1_V1} model, which we found by experimentation to match well -to the observed profiles generated by PS1. -Figure~\ref{fig:radial.profiles} shows example radial profiles for -moderately bright stars in fairly good (0.9 arcsec) and poor (2.2 -arcsec) seeing. Using a fixed power-law exponent results in somewhat -faster profile fitting compared to the variable power-law exponent -model. - -% moffat : 1969A&A.....3..455M -% buonanno : 1983A&AS...51...83B - -\begin{figure}[htbp] - \begin{center} - \includegraphics[width=\hsize]{{pics/radial.profiles}.\plotext} - \caption{\label{fig:radial.profiles} Radial profiles of stellar images from PS1. These two - profiles illustrate the radial trend of the PS1 PSFs for a star - with FWHM 0.9 arcsec (red) and 2.2 arcsec (blue). The black line - shows the PSF model with radial trend of the form $(1 + \kappa r^2 + r^{3.33})^{-1}$.} - \end{center} -\end{figure} - -\subsection{Radial Profiles} +\subsection{Extended Source Analysis} + +\note{intro paragraph: After the initial, fast analysis of the image + relying primarily on the PSF model, a complete analysis of the + extended source properties may be performed. For PS1 processing, + this step is the nightly (PV0) analysis of individual exposures and + only performed for the stacks. } + +\subsubsection{Radial Profiles} Galaxies with regular profiles, such as elliptical galaxies and @@ -1391,5 +1376,5 @@ % \note{these profiles are not saved in PSPS} -\subsection{Petrosian Radii and Magnitudes} +\subsubsection{Petrosian Radii and Magnitudes} \cite{1976ApJ...209L...1P} defined an adaptive aperture based on a @@ -1397,5 +1382,5 @@ aperture which can be determined for galaxies without significant biases as a function of distance. Since surface brightness in a -resolved object is conserved, using a ratio of surface brightness to +resolved source is conserved, using a ratio of surface brightness to define a spatial scale results in a spatial scale which is constant regardless of galaxy distance. @@ -1436,10 +1421,10 @@ Petrosian flux is contained. -\subsection{Radial Profile Wings} +\subsubsection{Radial Profile Wings} We attempt to measure the radial profile of sources in order to find -the radius at which the flux of the object is matches the sky. In +the radius at which the flux of the source is matches the sky. In this analysis, a series of up to 25 radial bins with power-law spacing -are defined and the flux of the object in each annulus is measured. +are defined and the flux of the source in each annulus is measured. The ``sky radius'' is defined to be the radius at which the (robust median) flux in the annulus is within 1 $\sigma$ of the local sky @@ -1451,5 +1436,6 @@ calculation of the kron magnitude. -\subsection{Kron Magnitudes} +\subsubsection{Kron Magnitudes} +\label{sec:kron.mags} Preliminary Kron radius and flux values \citep{1980ApJS...43..305K} @@ -1488,5 +1474,5 @@ opposites sides of the central pixel are considered together. The geometric mean of the two fluxes is used to replace the flux values. -If the object has 180\degree\ symmetry, this operation has no impact. +If the source has 180\degree\ symmetry, this operation has no impact. However, if one of the two pixels is unusually high, the value will be surpressed by the matched pixel on the other side. This trick has the @@ -1494,5 +1480,5 @@ neighbors. -\subsection{Convolved Galaxy Model Fits} +\subsubsection{Convolved Galaxy Model Fits} In the galaxy model fittting stage, sources which meet certain @@ -1526,5 +1512,5 @@ in the analysis. This restriction limited the total time spent on the galaxy modeling analysis at the expense of galaxy photometry in the -plane (though Kron photometry is available for those objects). The +plane (though Kron photometry is available for those sources). The Galactic Plane region was defined by $|b| > b_{\rm min}$ where $b_{\rm min} = b_0 + r_b e^{\frac{-l^2}{2 \sigma_b^2}}$. For the PV3 @@ -1662,5 +1648,5 @@ % DOI: https://doi.org/10.1071/AS05001 -\subsection{Convolved Radial Aperture Photometry} +\subsubsection{Convolved Radial Aperture Photometry} For some science goals, a well-measured color of a galaxy is more @@ -1676,9 +1662,9 @@ 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 +of the apertures for each source 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 determined based on the +for all 5 filters. In this analysis, the best model for each source +is subtracted from the image pixels for all sources excluding the +source in consideration. The 'best model' is determined based on the minimum $\chi^2$ value for the model fits. @@ -1689,5 +1675,5 @@ 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 +best source models are subtracted from the image for sources not being measured. This subtraction includes the convolution to smooth the model to the effective FWHM of the convolved image. The entire @@ -1695,4 +1681,192 @@ % \note{is the first convolution done with the Alard-Lupton technique?} + +\subsection{Aperture Correction} + +The important concept here is that an analytical model will always +fail to describe the flux of the sources at some level. In the end, +all astronomical photometry is in some sense a relative measurement +between two images. Whether the goal is calibration of a science +image taken at one location to a standard star image at another +location, or the goal is simply the repetitive photometry of the same +star at the same location in the image, it is always necessary to +compare the photometry between two images. If this measurement is to +be consistent, then the measurement must represent the flux of the +stars in the same way regardless of the conditions under which the +images were taken, at least within some range of normal image +conditions. So, for example, two images with different image quality, +or with different tracking and focus errors, will have different PSF +models. Since an analytical model will always fail to represent the +flux of the star at some level, the measured flux of the same source +in the two images will be different (even assuming all other +atmospheric and instrumental effects have been corrected). The +amplitude of the error will by determined by how inconsistently the +models represent the actual source flux. For example, if the first +image PSF model flux is consistently 10\% too low and the second is 5\% +too high, then the comparison between the two images will be in error +by 15\%. + +Aperture photometry avoids these problems, by trading for other +difficulties. In aperture photometry, if a large enough aperture is +chosen, the amount of flux which is lost will be a small fraction of +the total source flux. Even more importantly, as the image conditions +change, the amount lost will change by an even smaller fraction, at +least for a large aperture. This can be seen by the fact that the +dominant variations in the image quality are in the focus, tracking +and seeing. All of these errors initially affect the cores of the +stellar images, rather than the wide wings. The wide wings are +largely dominated by scattering in the optics and scattering in the +atmosphere. The amplitude and distribution of these two scattering +functions do not change significantly or quickly for a single +telescope and site. + +The difficulty for aperture photometry is the need to make an accurate +measurement of the local background for each source. As the aperture +grows, errors in the measurement of the sky flux start to become +dominant. If the aperture is too small, then variation in the image +quality are dominant. The brighter is the source, the smaller is the +error introduced by the large size of the aperture. However, the +number of very bright stars is limited in any image, and of course the +brighter stars are more likely to suffer from non-linearity or +saturation. PSPhot measures the aperture correction ({\em ApResid}) +for every PSF candidate source and applies this correction to the PSF +model photometry. + +% How important is this effect? Consider a typical bright source 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 source +% 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 +%%% candidate source, then calculates the trend of this correction as a +%%% function of the magnitude. This trend is fitted with a line. The +%%% resulting function can be used to determine the effective aperture +%%% correction for an infinite flux source and the average bias inherent +%%% in the sky measurement for the image. The scatter of the +%%% PSF-candidate source measurements about this trend is a measure of how +%%% well we can measure photometry from the image by applying the specific +%%% PSF model. The slope of this trend is a measure of the bias in the +%%% local sky measurment for each source. In principal, the measured sky +%%% levels could be modified by this bias. More generally, the measured +%%% bias in a collection of images could be used to improve the model +%%% fitting or sky fitting portion of the software the remove the bias +%%% term. + +PSPhot allows a collection of PSF model functions to be tried on all +PSF candidate sources. For each model test, the above corrected +ApResid scatter is measured. The PSF model function with the smallest +value for the ApResid scatter is then used by PSPhot as the best PSF +model for this image. The number of models to be tested is specified +by the configuration keyword \code{PSF_MODEL_N}. The configuration +variables \code{PSF_MODEL_0}, \code{PSF_MODEL_1}, through +\code{PSF_MODEL_N - 1} specify the names of the models which should be +tested. + +Several likely PSF model classes are available within \code{psphot}: +\begin{itemize} +\item Gaussian : $f = I_0 e^{-z}$ +\item Pseudo-Gaussian : $f = I_0 (1 + z + \frac{1}{2} z^2 + \frac{1}{6} z^3)^{-1}$ \code{[PGAUSS]} +\item Variable Power-Law : $f = I_0 (1 + z + z^{\alpha})^{-1}$ \code{[RGAUSS]} +\item Steep Power-Law : $f = I_0 (1 + \kappa z + z^{2.25})^{-1}$ \code{[QGAUSS]} +\item PS1 Power-Law : $f = I_0 (1 + \kappa z + z^{1.67})^{-1}$ \code{[PS1_V1]} +\end{itemize} +where $z \propto r^2$ ($z = \frac{x^2}{2\sigma_x^2} + +\frac{y^2}{2\sigma_y^2} + \sigma_{\rm xy} x y $). The Pseudo-Gaussian +is a Taylor expansion of the Gaussian and is used by Dophot +\citep{1993PASP..105.1342S}. The latter profiles are similar to the +Moffat profile form \citep{1969AA.....3..455M,1983AA...126..278B}, +with small differences. For the PS1 GPC1 analysis, we used the +\code{PS1_V1} model, which we found by experimentation to match well +to the observed profiles generated by PS1. +Figure~\ref{fig:radial.profiles} shows example radial profiles for +moderately bright stars in fairly good (0.9 arcsec) and poor (2.2 +arcsec) seeing. Using a fixed power-law exponent results in somewhat +faster profile fitting compared to the variable power-law exponent +model. + +% moffat : 1969A&A.....3..455M +% buonanno : 1983A&AS...51...83B + +\begin{figure}[htbp] + \begin{center} + \includegraphics[width=\hsize]{{pics/radial.profiles}.\plotext} + \caption{\label{fig:radial.profiles} Radial profiles of stellar images from PS1. These two + profiles illustrate the radial trend of the PS1 PSFs for a star + with FWHM 0.9 arcsec (red) and 2.2 arcsec (blue). The black line + shows the PSF model with radial trend of the form $(1 + \kappa r^2 + r^{3.33})^{-1}$.} + \end{center} +\end{figure} + +\subsection{Output Formats} + +\section{Forced Photometry Modes} + +\subsection{Forced Photometry : PSFs} + +\subsection{Forced Photometry : galaxies} + +\section{Difference Image Photometry} + +The variance map for a difference image must be generated from the two +images used to construct the difference. Otherwise, the low sky level +will automatically result in inconsistent interpretation of the variance. + +For a difference image, both positive and negative sources will be +present. The basic peak detection algorithm will only trigger for the +positive sources. One solution is to simply apply PSPhot to both the +difference image and its negative value. \note{do we want to code in +an automatic switch to get both positive and negative excursions in +the single pass?}. + +In the case of a difference image, the PSF model construction stage +will probably fail for lack of valid sources. It is better in these +cases to provide PSF model from some other source. For example, the +two images which are combined to generate the difference image +represent the PSF. Presumably, one or both have been convolved with a +PSF-matching kernel. The images which result from the convolution +should be used to measure the PSF model. + +The source classification scheme defaults to the galaxy models for +sources which are not well represented by the PSF model. In a +properly-constructed difference image, galaxies are unlikely to remain +behind as significant sources. Most real sources in the difference +image will be PSF-like and will consist of photometrically variable +sources (flare stars, supernovae, etc) or astrometrically variable +sources (high-proper motion stars or solar-system bodies). There are +three likely classes of sources which will not be well represented by +the PSF model. 1) Fast-moving solar-system objects will appear as +short streaks. For example, a fast solar system object would have an +apparent rate of 0.5 degrees per hour, translating to 15 arcseconds in +a 30 second exposure. Even a main belt asteroid at roughly 1 AU would +have reflect motion of approximately 1 degree per day, equivalent to +1.25 arcsec in a 30 second exposure, and could be noticeably smeared +and non-PSF-like. A trailed-star model can be used to characterize +these types of sourcess. 2) Small offset stars, either due to +atmospheric / color effects or modest proper motion will appear as PSF +dipoles in the difference images. The positive and the negative +images will have stellar profiles, but they will be significantly +offset and will not subtract well. The two components may not have +the same amplitude. A PSF-dipole model can be used to fit these types +of sources, with free parameters of the two centroids and the two +fluxes. 3) Comets will appear in the difference images as a non-PSF +sources. Their 2-D structure includes both the flux from the coma +(with a typical power-law profile) and flux from the tail (with a more +complex flux distribution). A comet flux model can be used to +characterize these sources in difference images. A major difficulty +in applying these three types of models is in making a robust test of +which model should be used. This problem is akin to the issue of +selecting and distinguishing between multiple galaxy models, as +discussed in the section on Galaxy models. + +\section{Examples and Tests} \acknowledgments @@ -1720,8 +1894,4 @@ \end{document} -\subsection{Forced Photometry : PSFs} - -\subsection{Forced Photometry : galaxies} - \subsection{Output Options} @@ -1729,65 +1899,4 @@ % \note{need to discuss failings and holes} - -\section{Alternative Scenarios} - -\subsection{Trailed Sources} - -\subsection{Difference Images} - -The variance map for a difference image must be generated from the two -images used to construct the difference. Otherwise, the low sky level -will automatically result in inconsistent interpretation of the variance. - -For a difference image, both positive and negative objects will be -present. The basic peak detection algorithm will only trigger for the -positive sources. One solution is to simply apply PSPhot to both the -difference image and its negative value. \note{do we want to code in -an automatic switch to get both positive and negative excursions in -the single pass?}. - -In the case of a difference image, the PSF model construction stage -will probably fail for lack of valid sources. It is better in these -cases to provide PSF model from some other source. For example, the -two images which are combined to generate the difference image -represent the PSF. Presumably, one or both have been convolved with a -PSF-matching kernel. The images which result from the convolution -should be used to measure the PSF model. - -The object classification scheme defaults to the galaxy models for -objects which are not well represented by the PSF model. In a -properly-constructed difference image, galaxies are unlikely to remain -behind as significant sources. Most real objects in the difference -image will be PSF-like and will consist of photometrically variable -objects (flare stars, supernovae, etc) or astrometrically variable -objects (high-proper motion stars or solar-system objects). There are -three likely classes of objects which will not be well represented by -the PSF model. 1) Fast-moving solar-system objects will appear as -short streaks. For example, a fast solar system object would have an -apparent rate of 0.5 degrees per hour, translating to 15 arcseconds in -a 30 second exposure. Even a main belt asteroid at roughly 1 AU would -have reflect motion of approximately 1 degree per day, equivalent to -1.25 arcsec in a 30 second exposure, and could be noticeably smeared -and non-PSF-like. A trailed-star model can be used to characterize -these types of objects. 2) Small offset stars, either due to -atmospheric / color effects or modest proper motion will appear as PSF -dipoles in the difference images. The positive and the negative -images will have stellar profiles, but they will be significantly -offset and will not subtract well. The two components may not have -the same amplitude. A PSF-dipole model can be used to fit these types -of objects, with free parameters of the two centroids and the two -fluxes. 3) Comets will appear in the difference images as a non-PSF -objects. Their 2-D structure includes both the flux from the coma -(with a typical power-law profile) and flux from the tail (with a more -complex flux distribution). A comet flux model can be used to -characterize these objects in difference images. A major difficulty -in applying these three types of models is in making a robust test of -which model should be used. This problem is akin to the issue of -selecting and distinguishing between multiple galaxy models, as -discussed in the section on Galaxy models. - -\subsection{Input \& Output Data Formats} - -\section{Sample Tests} \begin{verbatim} @@ -1827,5 +1936,5 @@ * authors * PSF residual map -* section 3.5.3 Model applied to detected objects needs to be reviewed +* section 3.5.3 Model applied to detected sources needs to be reviewed * read for english & phrasing