IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 22, 2016, 9:06:08 AM (10 years ago)
Author:
eugene
Message:

updates to algorithm descriptions

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/doc/release.2015/ps1.analysis/analysis.tex

    r39819 r39820  
    9696\section{INTRODUCTION}\label{sec:intro}
    9797
    98 \note{more PS1 background}
    99 
    10098The Pan-STARRS Image Processing Pipeline is responsible for the basic
    10199analysis of images from the Pan-STARRS telescopes Gigapixel Camera.
    … …  
    115113
    116114An additional constraint on the Pan-STARRS system comes from the high
    117 data rate.  PS1 produces typically $\sim 700$ GB per night of imaging
    118 data.  These images range from high galactic latitudes to the Galactic
    119 Bulge, so large numbers of measurable stars can be expected in much of
    120 the data.  The combination of the high precision goals of the
    121 astrometric and photometric measurements and the high data rate (and a
    122 finite computing budget) mean that the process of detecting,
    123 classifying, and measuring the astronomical objects in the image data
    124 stream in a timely fashion are a significant challenge.
     115data rate.  PS1 produces typically $\sim 500$ exposures per night,
     116corresponding to $\sim 750$ billion pixels of imaging data.  The
     117images range from high galactic latitudes to the Galactic bulge, so
     118large numbers of measurable stars can be expected in much of the data.
     119The combination of the high precision goals of the astrometric and
     120photometric measurements and the high data rate (and a finite
     121computing budget) mean that the process of detecting, classifying, and
     122measuring the astronomical objects in the image data stream in a
     123timely fashion are a significant challenge.
    125124
    126125In order to achieve these ambitious goals, the object detection,
    … …  
    154153  automated fashion, does it handle 2D variations well? (P. Stetson)
    155154
    156 \item Sextractor : pure aperture measurement with rudimentary
    157   object subtraction.  pro: fast, widely used, easy to automate.  con:
    158   poor object separation in crowded regions, PSF-modeling is only
    159   beta (psfex), what models are available? (E. Bertin)
     155\item Sextractor : pure aperture measurement with rudimentary object
     156  subtraction.  pro: fast, widely used, easy to automate.  con: poor
     157  object separation in crowded regions, PSF-modeling was only in beta,
     158  not widely used at the time. (E. Bertin)
    160159
    161160\item apphot : IRAF-based aperture photometry.  pro: widely used.
    … …  
    172171\end{itemize}
    173172
    174 \note{re-phrase this:} The Pan-STARRS IPP team decided that none of
    175 the existing packages met all of their needs, particularly given the
    176 very challenging goals of the project.  We decided to redesign the
    177 photometry analysis from scratch, using the lessons learned from the
    178 existing photometry systems.  In the process, the object analysis
    179 software would be written using the data analysis C-code library
    180 written for the IPP, \code{psLib}, and the components of the
    181 photometry code would be integrated into the IPP's mid-level astronomy
    182 data analysis toolkit called \code{psModules}.  The result is
    183 'PSPhot', which can be used either as a stand-alone C program, or as
    184 callable set of functions.
    185 
    186 \note{discuss the psphot program varients}
    187 
    188 \begin{verbatim}
    189 Other Varients:
    190 * psphotStack -- 5 filter simultaneous fitting
    191 * psphotFullForce
    192 \end{verbatim}
     173When the IPP development was starting, the existing photometry
     174packages either did not meet the level of accuracy required or were
     175required too much human intervention to be considered for the needs of
     176PS1.  In the case of the SDSS Photo tool, the software was judged to
     177be too tightly integrated to the architecture of SDSS to be easily
     178re-integrated into the Pan-STARRS pipelin.  A new photometry analysis
     179package was developed using lessons learned from the existing
     180photometry systems.  In the process, the object analysis software was
     181written using the data analysis C-code library written for the IPP,
     182\code{psLib}.  Components of the photometry code were integrated into
     183the IPP's mid-level astronomy data analysis toolkit called
     184\code{psModules}.  The result software, 'PSPhot', can be used either
     185as a stand-alone C program, or as a set of library functions which may
     186be integrated into other programs
     187
     188The main version of PSPhot is a stand-alone program which is run on a
     189single image or a group of related images representing the data read
     190from a camera in a single exposure.  The images are expected to have
     191already been detrended so that pixel values are linearly related to
     192the flux.  The gain may be specific by the configuration system, or a
     193variance image may be supplied.  A mask may also be supplied to mark
     194good, bad, and suspect pixels.  Several variants of psphot have also
     195been used in the PS1 PV3 analysis. 
     196
     197The version called PSPhotStack accepts a set of images, each
     198representing the same patch of sky in a different filter, nominally
     199the full $grizy$ filter set for the analysis of the PS1 PV3 stack
     200images, though where insufficient data were available in a given
     201filter, a subset of these filters was processed as a group.  As
     202discussed in detail below, the PSPhotStack analysis includes the
     203capability of measuring forced PSF photometry in some filter images
     204based on the position of sources detected in the other filters.  It
     205also include an option to convolve the set of images to a single,
     206common PSF size across the filters for the purpose of fixed aperture
     207photometry.
     208
     209A second version of PSPhot used in the PV3 analysis is called
     210PSPhotFullForce.  In this version, a set of image all representing the
     211same pixels are processed together, with the positions of sources to
     212be analysed loaded from a supplied file.  In this version the
     213analysis, sources are not discovered -- only the supplied sources are
     214considered.  PSF models are determined for each exposure and the
     215forced PSF photometry is measured for all sources.  A subset of
     216sources may also be used to measure forced galaxy shape parameters.
     217As described below, a grid of galaxy models are fitted based on the
     218supplied guess model. 
    193219
    194220\section{PSPhot Design Goals}
    … …  
    224250  level must be reached for images with 250 mas pixels, implying
    225251  PSPhot must introduce measurement errors less than 1/50th of a
    226   pixel. \note{the choice of F32 parameters places a numerical limit
    227   of 1e-7 on the accuracy of a pixel relative to the size of a chip
    228   (since a single data value is used for X or Y).  For the $4800^2$
    229   GPC chips, this yields a limit of about 0.25 milliarcsecond.}
     252  pixel. The choice of 32 bit floating point data values for the
     253  source centroids places a numerical limit of 1e-7 on the accuracy of
     254  a pixel relative to the size of a chip (since a single data value is
     255  used for X or Y).  For the $4800^2$ GPC chips, this yields a limit
     256  of about 0.25 milliarcsecond.
    230257\end{itemize}
    231258
    … …  
    244271
    245272\item {\bf Flexible non-PSF models} PSPhot must be able to represent
    246   PSF-like objects as well as non-PSF sources.  It must be easy to add
    247   new object models as interesting representations of sources are
    248   invented.
     273  PSF-like objects as well as non-PSF sources (e.g., galaxies).  It
     274  must be easy to add new object models as interesting representations
     275  of sources are invented.
    249276
    250277\item {\bf Clean code base} PSPhot should incorporate a high-degree of
    … …  
    256283  provide the user with methods for assessing the different PSF models.
    257284
    258 \item {\bf Careful aperture corrections} PSPhot must carefully measure
    259   and correct for the photometric and astrometric trends introduced by
    260   using analytical PSF models.
     285\item {\bf Careful systematic corrections} PSPhot must carefully
     286  measure and correct for the photometric and astrometric trends
     287  introduced by using analytical PSF models.
    261288
    262289\item {\bf User Configurable} PSPhot should allow users to change the
    … …  
    276303
    277304\item {\bf Initial object detection} Smooth, find peaks, measure basic
    278   properties
     305  properties.
    279306
    280307\item {\bf PSF determination} Select PSF candidates, perform model
    … …  
    288315  properties (aperture or PSF)
    289316
     317\item {\bf Extended Source Analysis} Detailed measurements relevant to
     318  galaxies and/or other extended (non-PSF) sources.
     319
    290320\item {\bf Aperture corrections} Measure the curve-of-growth, spatial
    291321  aperture variations, and background-error corrections. 
    … …  
    296326
    297327PSPhot is highly configurable.  Users may choose via the configuration
    298 system which of the above analyses are performed.  This may be useful
    299 for testing, but may also allow for specialized use cases.  For
    300 example, the PSF model may already be available from external
    301 information, in which case the PSF modeling stage can be skipped.
     328system which of the above analyses are performed.  This is useful for
     329testing, but also allows for specialized use cases.  For example, the
     330PSF model may already be available from external information, in which
     331case the PSF modeling stage can be skipped.
    302332
    303333\subsection{Image Preparation}
    … …  
    343373circumstance, while a pixel in which persistence ghosts have been
    344374subtracted might be useful for detection or even analysis of brighter
    345 sources.  \note{can I identify which functions respect which sets of masks}
     375sources.  Table~\ref{tab:mask_values} lists the 16 bit values used for
     376PS1 mask images, along with their description (see \note{Waters et
     377  al. paper} for additional information).
     378
     379\begin{table}
     380\caption{\label{tab:mask_values} PSPhot / GPC1 Mask Image Pixel Values}\vspace{-0.5cm}
     381\begin{center}
     382\begin{tabular}{lcl}
     383\hline
     384\hline
     385{\bf Mask Name} & {\bf Mask Value} & {\bf Description} \\
     386\hline
     387  DETECTOR & 0x0001 & A detector defect is present. \\
     388  FLAT     & 0x0002 & The flat field model does not calibrate the pixel reliably. \\
     389  DARK     & 0x0004 & The dark model does not calibrate the pixel reliably. \\
     390  BLANK    & 0x0008 & The pixel does not contain valid data. \\
     391  CTE      & 0x0010 & The pixel has poor charge transfer efficiency. \\
     392  SAT      & 0x0020 & The pixel is saturated. \\
     393  LOW      & 0x0040 & The pixel has a lower value than expected. \\
     394  SUSPECT  & 0x0080 & The pixel is suspected of being bad. \\
     395  BURNTOOL & 0x0080 & The pixel contain an burntool repaired streak. \\
     396  CR       & 0x0100 & A cosmic ray is present. \\
     397  SPIKE    & 0x0200 & A diffraction spike is present. \\
     398  GHOST    & 0x0400 & An optical ghost is present. \\
     399  STREAK   & 0x0800 & A streak is present. \\
     400  STARCORE & 0x1000 & A bright star core is present. \\
     401  CONV.BAD & 0x2000 & The pixel is bad after convolution with a bad pixel. \\
     402  CONV.POOR& 0x4000 & The pixel is poor after convolution with a bad pixel. \\
     403  MARK     & 0x8000 & An internal flag for temporarily marking a pixel. \\
     404\hline
     405\end{tabular}
     406\end{center}
     407\end{table}
    346408
    347409The variance image, if not supplied is constructed by default from the
    … …  
    367429Since a typical smoothing or warping operation may introduce
    368430correlation between 25 - 100 neighboring pixels, the size of such a
    369 covariance image is prohibitive.  In practice, however, there are two
    370 extreme cases which generally are relevant.  \note{talk about the
    371   covar matrix for a PSF}
    372 
    373 \subsection{Background (Sky) Model}
     431covariance image is prohibitive. 
     432
     433Before sources are detected in the image, a model of the background is
     434subtracted.  The image is divided into a grid of background points
     435with a spacing of 400 pixels.  Superpixels of size $800\times 800$
     436pixels are used to measure the local background for each background
     437grid point, thus over-sampling the background spatial variations by a
     438factor of 2.  In the interest of speed, 10,000 randomly selected
     439{\em unmasked} pixels in these regions are sampled to determine the
     440background.  \note{flesh out the details here}.  Bilinear
     441interpolation is used to generate a full-resolution image from the grid of
     442background points, and this image is then subtracted from the science
     443image.  The background image and the background standard deviation
     444image are kept in memory from which the values of \code{SKY} and
     445\code{SKY\_SIGMA} are calculated for each object in the output catalog.
    374446
    375447\subsection{Initial Object Detection}
     448
     449\subsubsection{Peak Detection}
     450\label{sec:peaks}
    376451
    377452The objects are initially detected by finding the location of local
    378453peaks in the image.  The flux and variance images are smoothed with a
    379 small circularly symmetric kernel using a two-pass 1D Gaussian
    380 (\note{KEYWORD?}).  The smoothed flux and variance images are combined
    381 to generate a significance image in signal-to-noise units
    382 \note{including correction for the covariance, if known}. At this
    383 stage, the goal is only to detect the brighter sources, above a user
    384 defined S/N limit (configuration keyword: \code{PEAK\_NSIGMA}).  The
    385 detection efficiency for the brighter sources is not strongly
    386 dependent on the form of this smoothing function.
     454small circularly symmetric kernel using a two-pass 1D Gaussian.  The
     455smoothed flux and variance images are combined to generate a
     456significance image in signal-to-noise units, including correction for
     457the covariance, if known. At this stage, the goal is only to detect
     458the brighter sources, above a user defined S/N limit (configuration
     459keyword: \code{PEAKS\_NSIGMA\_LIMIT}).  A maximum of
     460\code{PEAKS\_NMAX} are found at this stage.  The detection efficiency
     461for the brighter sources is not strongly dependent on the form of this
     462smoothing function.
    387463
    388464The local peaks in the smoothed image are found by first detecting
    … …  
    397473the maximum $X$ and $Y$ corners of the region.
    398474
    399 \subsection{Footprints}
    400 
    401 \note{need to describe the process of generating the source footprints
    402   and then culling the insignificant peaks}
    403 
    404 \subsubsection{Moments and related}
    405 
    406 \note{disucss the Kron mags}
    407 
    408 \note{this section is wrong: we no longer use S/N clipping, but a
    409   Gaussian window function, chosen based on the measured moment}
    410 
    411 Once a collection of peaks have been identified, basic properties of
    412 the objects are measured.  First, the local sky flux is measured
    413 within a square annulus with user-defined dimensions
    414 (\code{INNER\_RADIUS} and \code{OUTER\_RADIUS}), using the sample
    415 median.  This local background value is then used to calculate the
    416 object first and second moments within a small user-defined aperture
    417 (\code{MOMENT\_RADIUS}).  The first-order moments are a good
    418 representation of the object position, while the second-order moments
    419 are a measure of the object shape.  The second-order moments are
    420 somewhat sensitive to the size of the aperture and the accuracy of the
    421 background measurement.  The moment calculation is only performed
    422 using pixels which exceed a S/N of 1.  If, in the process of
    423 calculating the source moments, the S/N limits reject all but \note{3}
    424 or fewer of the source pixels, the peak is identified as being
    425 suspect, and is not used for further analysis.  If the measured
    426 centroid coordinates differ from the peak coordinates be a large
    427 amount (\code{MOMENT\_RADIUS}), then the peak is again identified as
    428 being of poor quality and is rejected.  In both of these cases, it is
    429 likely that the `peak' was identified in a region of flat flux
    430 distribution or many saturated or edge pixels.
    431 
    432 \subsubsection{Determination of the Peak Coordinates and Errors}
    433 
    434475We use the 9 pixels which include the source peak to fit for the
    435476position and position errors.  We model the peak of the sources as a
    … …  
    448489of only 0 or 1, we can greatly simplify the chi-square equation to a
    449490square matrix equation with the following values:
    450 
    451 %% fix this:
    452 \begin{verbatim}
    453 | 9 0 0 0 6 6 | C_00 | = \sum F_{i,j}
    454 | 0 6 0 0 0 0 | C_10 | = \sum F_{i,j} x
    455 | 0 0 6 0 0 0 | C_01 | = \sum F_{i,j} y
    456 | 0 0 0 6 0 0 | C_11 | = \sum F_{i,j} x y
    457 | 6 0 0 0 6 4 | C_20 | = \sum F_{i,j} x^2
    458 | 6 0 0 0 4 6 | C_02 | = \sum F_{i,j} y^2
    459 \end{verbatim}
    460 
    461 The inverse of the 3x3 matrix terms for $C_{00}$, $C_{20}$, and $C_{02}$ is:
    462 \begin{verbatim}
    463 | +5/9 -1/3 -1/3 |
    464 | -1/3 +1/2    0 |
    465 | -1/3    0 +1/2 |
    466 \end{verbatim}
    467 
    468 The location of the peak is determined from the minimum of the
     491\[
     492\left( \begin{array}{cccccc}
     4939 & 0 & 0 & 0 & 6 & 6 \\
     4940 & 6 & 0 & 0 & 0 & 0 \\
     4950 & 0 & 6 & 0 & 0 & 0 \\
     4960 & 0 & 0 & 6 & 0 & 0 \\
     4976 & 0 & 0 & 0 & 6 & 4 \\
     4986 & 0 & 0 & 0 & 4 & 6 \\
     499\end{array} \right)
     500\left( \begin{array}{c}
     501C_{00}\\
     502C_{10}\\
     503C_{01}\\
     504C_{11}\\
     505C_{20}\\
     506C_{02}\\
     507\end{array} \right)
     508=
     509\left( \begin{array}{c}
     510\sum F_{i,j}     \\
     511\sum F_{i,j} x   \\
     512\sum F_{i,j} y   \\
     513\sum F_{i,j} x y \\
     514\sum F_{i,j} x^2 \\
     515\sum F_{i,j} y^2 \\
     516\end{array} \right)
     517\]
     518
     519Inverting the 3x3 matrix terms for $C_{00}$, $C_{20}$, and $C_{02}$,
     520the location of the peak is determined from the minimum of the
    469521bi-quadratic function above, and is given by:
    470 
    471522\begin{eqnarray}
    472 Det    & = & 4 C_{20} C_{02} - C_{11}^2 \\
    473 x_{min} & = & (C_{11} C_{01} - 2 C_{02} C_{10}) / Det \\
    474 y_{min} & = & (C_{11} C_{10} - 2 C_{20} C_{01}) / Det \\
     523x_{min} & = & (C_{11} C_{01} - 2 C_{02} C_{10}) D^{-1} \\
     524y_{min} & = & (C_{11} C_{10} - 2 C_{20} C_{01}) D^{-1} \\
     525D      & = & 4 C_{20} C_{02} - C_{11}^2
    475526\end{eqnarray}
    476527
    477 \note{error on the peak position}
     528\subsubsection{Footprints}
     529
     530The peaks detected in the image may correspond to real sources, but
     531they may also correspond to noise fluctuations, especially in the
     532wings of bright stars.  PSPhot attempts to identify peaks which may be
     533formally significant, but are not locally significant.  It first
     534generates a set of ``footprints'', contiguous collections of pixels in
     535the smoothed significance image above the detection threshold.  These
     536regions are grown by a small amount to avoid errors on rough edges --
     537an image of the footprints is convolved with a disk of radius 3
     538pixels.  Peaks are assigned to the footprints in which they are
     539contained (note by definition all peaks must be located in a
     540footprint). 
     541
     542For any peak which is not the brightest peak in that footprint it is
     543possible to reach the brightest peak by following the highest valued
     544pixels between the two peaks.  The lowest pixel along this path is the
     545{\em key col} for this peak (as used in topographic descriptions of a
     546mountain).  If the key col for a given peak is less than
     547\code{FOOTPRINT\_CULL\_NSIGMA\_DELTA} (4.0) sigmas below the peak of
     548interest, the peak is considered to be {\em locally insignificant} and
     549removed from the list of possible detections.  In the vicinity of a
     550saturated star, the rule is somewhat more agressive as the flat-topped
     551or structured saturated top of a bright star may appear as multiple
     552peaks with highly significant cols between them.  However, this is an
     553artifact of the proximity to saturation.  In this regime, we require
     554the col to also be a fixed fraction (5\%) of the saturation below the
     555peak to avoid being marked as locally insignificant.
     556
     557\subsubsection{Centroid and higher-order Moments}
     558
     559Once a collection of peaks has been identified, a number of basic
     560properties of the objects related to the first and second moments are
     561measured.  Below, the second moments are used to select candidate
     562stellar sources to be used in modeling the PSF.
     563
     564In order to measure the moments, it is necessary to define an
     565appropriate aperture in which the moments are measured.  We also apply
     566a ``window function'', down-weighting the pixels by a Gaussian of size
     567$\sigma_W$ which is chosen to be large compared to the PSF size.  The
     568choice of the window function $\sigma_W$ and the aperture is an
     569iterative process: for a given value of $\sigma_W$, the PSF stars will
     570have a measured value of $\sigma$ which is smaller than the true value
     571due to the window function.  \note{generate examples to illustrate
     572  this}.
     573
     574To choose the value of $sigma_W$, we try values of (1, 2, 3, 4.5, 6,
     5759, 12, 18) pixels.  For each of these values, we then select candidate
     576PSF stars based on the distribution of the measured $\sigma_{x,x},
     577\sigma_{y,y}$ values.  For each test value of $\sigma_w$, determine
     578the ratio $f = \frac{\sigma_{x,x} + \sigma{y,y}}{2 \sigma_w}$, i.e.,
     579the ratio of the window size to the observed PSF size.  We interpolate
     580to find a value of $\sigma_W$ for which $f$ is expected to be 0.65.
     581\note{what is the expected ratio of $\sigma_x$ to the true value?}.
     582We call this value the \code{MOMENTS\_GAUSS\_SIGMA}.  We use an
     583aperture with a radius of \code{PSF\_MOMENTS\_RADIUS} = 4$\times$
     584\code{MOMENTS\_GAUSS\_SIGMA} to select the pixels for the measurement.
     585
     586Once \code{PSF\_MOMENTS\_SIGMA} has been determined, moments are
     587measured as defined below. 
     588
     589\begin{eqnarray}
     590x_0      & = & \frac{1}{S} \sum_i (f_i - s_i)x_i w_i \\
     591y_0      & = & \frac{1}{S} \sum_i (f_i - s_i)y_i w_i \\
     592M_{xx}   & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^2w_i \\
     593M_{xy}   & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)(y_i - y_0)w_i \\
     594M_{yy}   & = & \frac{1}{S} \sum_i (f_i - s_i)(y_i - y_0)^2w_i \\
     595M_{xxx}  & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^3w_i / r_i \\
     596M_{xxy}  & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^2(y_i - y_0)w_i / r_i \\
     597M_{xyy}  & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)(y_i - y_0)^2w_i / r_i \\
     598M_{yyy}  & = & \frac{1}{S} \sum_i (f_i - s_i)(y_i - y_0)^3w_i / r_i \\
     599M_{xxxx} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^4w_i / r^2_i \\
     600M_{xxxy} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^3(y_i - y_0)w_i / r^2_i \\
     601M_{xxyy} & = & \frac{1}{S} \sum_i (f_i - s_i)(x_i - x_0)^2(y_i - y_0)^2w_i / r^2_i \\
     602M_{xyyy} & = & \frac{1}{S} \sum_i (f_i - s_i)(y_i - y_0)(y_i - y_0)^3w_i / r^2_i \\
     603M_{yyyy} & = & \frac{1}{S} \sum_i (f_i - s_i)(y_i - y_0)^4w_i / r^2_i
     604\end{eqnarray}
     605where $f_i$ is the flux in a pixel; $s_i$ is the local sky value for
     606that pixel; $w_i$ is the value of the window function for the pixel;
     607$S = \sum_i (f_i - s_i) w_i$ is the window-weighted sum of the source
     608flux, used to re-normalize the moments; $r_i$ is the radius of a
     609pixel, $\sqrt{(x_i - x_0)^2 + (y_i - y_0)^2}$; The sum is performed
     610over all pixels in the aperture.  For the centroid calculation ($x_0,
     611y_0$), the peak coordinate (see~\ref{sec:peaks}) is used to define the
     612aperture and the window function; for higher order moments, the
     613centroid is used to center the window function.
     614
     615If the measured centroid coordinates ($x_0, y_0$) differs from the
     616peak coordinates be a large amount (\code{MOMENT\_RADIUS}), then the
     617peak is identified as being of poor quality and is rejected.  In
     618both of these cases, it is likely that the `peak' was identified in a
     619region of flat flux distribution or many saturated or edge pixels.
     620
     621In addition to the moments above, a preliminary Kron radius and flux
     622are also calculated at this stage.  In this analysis, the 1st and
     623half-radial moments are calculated:
     624\begin{eqnarray}
     625M_r & = & \frac{1}{S} \sum_i (f_i - s_i)r_i \\
     626M_h & = & \frac{1}{S} \sum_i (f_i - s_i)\sqrt{r_i}
     627\end{eqnarray}
     628Note that the window function is not applied in the calculation of
     629these moments.
     630
     631The Kron radius is defined the be 2.5$\times$ the first radial moment.
     632The Kron flux is the sum of (sky-subtracted) pixel fluxes within the
     633Kron radius.  We also calculate the flux in two related annular
     634apertures: the Kron inner flux is the sum of pixel values for the
     635annulus $R_1 < r < 2.5 R_1$, while the Kron outer flux is the sum of
     636pixel values for $2.5 R_1 < r < 4 R_1$.  The first radial moment is
     637limited at the low and high ends by $R_{\rm min} < M_r < R_{\rm max}$
     638where $R_{\rm min}$ is the first radial moment of the PSF stars, or
     6390.75$\times$ \code{MOMENTS\_GAUSS\_SIGMA} if that cannot be
     640determined.  $R_{\rm max}$ is set to \code{PSF\_MOMENTS\_RADIUS}, the
     641size of the moments aperture.
    478642
    479643\subsection{PSF Determination}
    … …  
    491655\begin{eqnarray}
    492656f(x,y) & = & I_o exp (-z) + S  \\
    493     R  & = & \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y -
    494     y_o)^2}{2\sigma_y^2} + \sigma_{\rm xy}(x - x_o)(y - y_o)
     657    z  & = & \frac{x^2}{2\sigma_x^2} + \frac{y^2}{2\sigma_y^2} + \sigma_{\rm xy} x y \\
     658    x  & = & x_{\rm ccd} - x_o \\
     659    y  & = & y_{\rm ccd} - y_o \\
    495660\end{eqnarray}
    496661The object model will have a variety of model parameters, in this case
    … …  
    503668The point-spread-function (PSF) of an image describes the shape of all
    504669unresolved objects in the image.  In a typical image, the shape of
    505 point sources is not well described by a single functional form;
    506 rather, the shape will vary as a function of position in the image.
    507 The PSF model therefore must describe the parameter variation as a
    508 function of the position of the object on the image.  Note that the
    509 object model consists of a certain number of parameters which are
    510 defined by the PSF model, and another set of parameters which are
    511 independent from object to object.  For the case of the elliptical
    512 Gaussian model, the PSF parameters would be the shape terms
    513 ($\sigma_x, \sigma_y, \sigma_{\rm xy}$) while the independent
    514 parameters would be the centroid, normalization and local sky values
    515 ($x_o, y_o, I_o, S$).  PSPhot uses a 2-D polynomial to specify the
    516 variation in the PSF parameters as a function of position in the
    517 image.  In the case of the elliptical Gaussian, this implies that the
    518 parameters are each a function of the object centroid coordinates:
     670point sources is not well described by a single function.  Instead,
     671the shape will vary as a function of position in the image.  The PSF
     672model therefore must describe the parameter variation as a function of
     673the position of the object on the image.  Note that the object model
     674consists of a certain number of parameters which are defined by the
     675PSF model, and another set of parameters which are independent from
     676object to object.  For the case of the elliptical Gaussian model, the
     677PSF parameters would be the shape terms ($\sigma_x, \sigma_y,
     678\sigma_{\rm xy}$) while the independent parameters would be the
     679centroid, normalization and local sky values ($x_o, y_o, I_o, S$).
     680PSPhot uses a 2-D polynomial to specify the variation in the PSF
     681parameters as a function of position in the image \note{or an
     682  interpolated map}.  In the case of the elliptical Gaussian, this
     683implies that the parameters are each a function of the object centroid
     684coordinates:
    519685\begin{eqnarray}
    520686\sigma_x    & = & f_1(x,y) \\
    … …  
    579745and the PSF model parameters.  It also defines a specific order for
    580746the 4 independent parameters. 
    581 
    582 \note{the code may also require that two of the PSF-like parameters
    583 represent the shape in some way}.
    584747
    585748\subsubsection{PSF Candidate Object Selection}
    … …  
    626789parameters are far from the minimization values.  PSPhot uses the
    627790first and second moments to make a good guess for the centroid and
    628 shape parameters for the PSF models.  \note{still true? In order to
    629   minimize the impact of close neighbors, the variance values used in
    630   the fit are enhanced by a fraction of the deviation of the
    631   particular pixel value from the model guess.}  Any objects which
    632 fail to converge in the fit are flagged as invalid.
    633 
    634 \note{does the variance enhancement introduce too much bias?}
    635 
    636 \note{discuss the convergence criteria, model parameter guesses}
     791shape parameters for the PSF models.  Any objects which fail to
     792converge in the fit are flagged as invalid.
    637793
    638794For the resulting collection of object model parameters, the
    639795PSF-dependent parameters of the models are all fitted as a function of
    640 position to a 2-D polynomial.  The order of this polynomial is (should
    641 be?) a user-defined parameter.  The fitting process for these
    642 polynomials is iterative, and rejects the $3-\sigma$ outliers in each
    643 of three passes.  This fitting technique results in a robust
    644 measurement of the variation of the PSF model parameters as a function
    645 of position without being excessively biased by individual objects
    646 which fail drastically.  Objects whose model parameters are rejected
    647 by this iterative fitting technique are also marked as invalid and
    648 ignored in the later PSF model fitting stages.
     796position to a 2-D polynomial.  The order of this polynomial is a
     797user-defined parameter.  The fitting process for these polynomials is
     798iterative, and rejects the $3-\sigma$ outliers in each of three
     799passes.  This fitting technique results in a robust measurement of the
     800variation of the PSF model parameters as a function of position
     801without being excessively biased by individual objects which fail
     802drastically.  Objects whose model parameters are rejected by this
     803iterative fitting technique are also marked as invalid and ignored in
     804the later PSF model fitting stages.
    649805
    650806All of the PSF-candidate objects are then re-fitted using the PSF
    … …  
    667823correction is judged to be the best model.
    668824
    669 \subsection{Very Bright Stars}
     825\subsection{Bright Source Analysis}
     826
     827\subsubsection{Very Bright Stars}
    670828\note{flesh out}
    671829
    … …  
    675833and subtract a radial profile modeled with a spline (?).
    676834
    677 \subsection{PSF vs CR vs Extended}
    678 
    679 \subsection{Bright Source Analysis}
    680 
    681 After a PSF model has been determined, PSPhot performs the analysis of
    682 the bright objects in the image.  In this stage, all of the objects
    683 with an estimated signal to noise (based on the moments analysis)
    684 greater than a user-set threshold are analysed and subtracted from the
    685 image.  An optional successive stage then finds fainter sources and
    686 measures them as well (see Faint Source Analysis,
    687 Section~\ref{faintsources}).  In the bright source analysis stage, two
    688 major varients are available.  In the primary version, all objects are
    689 examined (in decending order of brightness) and an appropriate models
    690 is determined for each object which is then subtracted; in the
    691 alternate version, the objects are examined (in decending order of
    692 brightness) and the PSF-like objects subtracted first, then the
    693 extended objects are analysed on a second pass.
    694 
    695835\subsubsection{Fast Ensemble PSF Fitting}
    696836
    … …  
    698838convenient to subtract off all of the sources, at least as well as
    699839possible at this stage.  We make the assumption that all sources are
    700 PSF-like.  We also assume their position can be taken as the peak of a
    701 2D quadratic function fitted to the peak pixel and its surrounding 8
    702 pixels.  A single linear fit is used to simultaneously measure all
    703 source fluxes.  Since the local sky has been subtracted, this
    704 measurement assumes the local sky is zero. 
    705 
     840PSF-like.  If the centroid of the source has been determined, we use
     841this value for its position; otherwise, we use the interpolated
     842position of the peak. A single linear fit is used to simultaneously
     843measure all source fluxes.  Since the local sky has been subtracted,
     844this measurement assumes the local sky is zero.  We can write a single
     845$\chi^2$ equation for this image:
    706846\[
    707 \chi^2 = \sum_{\rm pixels} (F_{x,y} - \sum_{\rm sources} A_i PSF[x,y])^2
     847\chi^2 = \sum_{\rm pixels} (F_{x,y} - \sum_{\rm sources} A_i P[x_0,y_0])^2
    708848\]
     849where $F_{x,y}$ is image flux for each pixel, $P[x_0,y_0]$ is the PSF
     850model realized at the position of source $i$, and $A_i$ is the
     851normalization for the source.
    709852
    710853Minimizing this equation with respect to each of the $A_i$ values
    … …  
    712855\[ M_{i,j} \bar{A_i} = \bar{F_j}\]
    713856where $\bar{A_i}$ is the vector of $A_i$ values, the elements of
    714 $M_{i,j}$ consist of the dot product of the unit-flux PSF for source
    715 $i$ and source $i$, and $\bar{F_j}$ is the dot product of the
    716 unit-flux PSF for source $i$ with the pixels corresponding to source
    717 $i$.  The dot products are calculated only using pixels within the
     857$M_{i,j}$ consist of the dot products of the unit-flux PSF for source
     858$i$ and source $j$, and $\bar{F_j}$ is the dot product of the
     859unit-flux PSF for source $j$ with the pixels corresponding to source
     860$j$.  The dot products are calculated only using pixels within the
    718861source apertures.  Since most sources have no overlap with most other
    719862sources, this matrix equation results in a very sparse, mostly
    720863diagonal square matrix.  The dimension is the number of sources,
    721 likely to be 1000s or 10,000s.  Such a matrix does not lend itself to
    722 direct inversion.  However, an interative solution quickly yields a
    723 result with sufficient accuracy.  In the iterative solution, a guess
    724 at the solution is made; the guess is multiplied by the matrix, and
    725 the result compared with the observed vector $\bar{F_j}$.  The
    726 difference is used to modify the initial guess. This proces is
    727 repeated several times to achieve a good convergence. 
     864likely to be 1000s or 10,000s.  Direct inversion of the matrix would
     865be computationally very slow.  However, an interative solution quickly
     866yields a result with sufficient accuracy.  In the iterative solution,
     867a guess at the solution $\bar{A}$ is made assuming $M_{i,j}$ is purely
     868diagonal; the guess is multiplied by $M_{i,j}$, and the result
     869compared with the observed vector $\bar{F_j}$.  The difference is used
     870to modify the initial guess.  This proces is repeated several times to
     871achieve a good convergence.
    728872
    729873Once a solution set for $A_i$ is found, all of the objects are
    … …  
    744888quality of the PSF model as a representation of the object shape.  To
    745889do this, it calculates the next step of the minimization {\em allowing
    746 the shape parameters to vary}.  This step, essentially the
     890  the shape parameters to vary}.  This step, essentially the
    747891Gauss-Newton minimization distance from the current local minimum,
    748892should be very small if the object is well represented by the PSF, but
    … …  
    752896elliptical Gaussian, these two parameters are $\sigma_x$ and
    753897$\sigma_y$.  For a generic model, the shape parameters may be defined
    754 differently, but the should always be two parameters which scale the
    755 object size in two dimensions (what about a polar-coordinate form?)
    756 Currently, PSPhot requires the two relevant shape parameters to be the
    757 first two dependent parameters in the list of model parameters (ie,
    758 parameters 4 \& 5).
     898differently, but there should always be two parameters which scale the
     899object size in two dimensions.  Currently, PSPhot requires the two
     900relevant shape parameters to be the first two dependent parameters in
     901the list of model parameters (ie, parameters 4 \& 5).
    759902
    760903The expected distribution of the variation of the two shape parameters
    761904will be a function of the signal-to-noise of the object in question
    762905and the value of the shape parameter itself.  The expected standard
    763 deviation on the shape parameter is, eg, $\sigma_x / {\rm SN}$.  If
     906deviation on the shape parameter is, eg, $\sigma_x / {\rm S/N}$.  If
    764907the object is well-represented by the PSF, then the shape parameter
    765908values should be close to their minimization value.  We can thus ask,
    … …  
    802945distribution, the remaining flux should be below 1 $\sigma$
    803946significance.  In practice the cores of bright stars are poorly
    804 represented and may have larger residual significance. \note{in future
    805 work, we may choose to enhance the variance to minimize detection of
    806 objects in the residuals of brighter objects}.
     947represented and may have larger residual significance.
     948
     949\note{I am not sure the above discussion is still (PV3) true.  To be reviewed.}
    807950
    808951\subsubsection{Blended Sources}
    … …  
    828971available non-PSF model or models.
    829972
    830 \note{better description of the acceptance criteria; the FLT model is
    831   tried before the DBL is accepted or rejected}.
     973\subsubsection{Source Size Assessment}
     974
     975After the PSF model has been fitted to all sources, and the Kron flux
     976has been measured for all sources, PSPhot uses these two measurements,
     977along with some additional pixel-level analysis, to determine the size class
     978of the object.  If the object is large compared to a PSF, it is
     979considered to be {\em extended} and will be
     980fitted with a galaxy model (or possibly another type of extended
     981source model in special cases).  If the object is small compared to a
     982PSF, it is considered to be a {\em cosmic ray} and masked.
     983
     984Extended sources are identified as those for which the Kron magnitude is
     985significantly brighter than the PSF magnitude when compared to a PSF
     986star.  The value $dMagKP = m_{\rm Kron} - m_{\rm PSF}$, the difference between the PSF
     987and Kron magnitudes, is calculated for each object.  The median of
     988$dMagKP$ is calculated for the PSF stars.  This median is subtracted
     989from $dMagKP$ for each star.  The result is divided by the quadrature
     990error of the PSF and Kron magnitudes and called \code{extNsigma}.  If
     991\code{extNsigma} is larger than \code{PSPHOT.EXT.NSIGMA.LIMIT} (3.0),
     992the object is considered to be extended.
     993
     994Cosmic Rays are identified by a combination of the Kron magnitude and
     995the second-moment width of the object in the minor axis direction.
     996The second-moment in the minor axis direction is calculated from
     997$M_{xx}, M_{xy}, M_{yy}$ as follows:
     998\[
     999M_{\rm minor} = \frac{1}{2}(M_{xx} + M_{yy}) - \frac{1}{2}\sqrt{(M_{xx} - M_{yy})^2 + 4 M_{xy}^2}
     1000\]
     1001If $M_{\rm minor} < 1.2$ pixels$^2$ and the instrumental Kron
     1002magnitude is $< -5.5$, then the object is identified as a cosmic ray
     1003and the associated pixels are masked.
    8321004
    8331005\subsubsection{Non-PSF Objects}
    … …  
    8731045\subsection{Faint Sources}
    8741046
    875 \note{this is not done : should use the ensemble PSF fitting to fit
    876   just the new significant peaks}
    877 
    8781047After a first pass through the image, in which the brighter sources
    8791048above a high threshold level have been detected, measured, and
    … …  
    8871056The objects which are measured in this faint-object stage are clearly
    8881057low significance detections.  A typical threshold for the bright
    889 object analysis would S/N of 5 - 10.  The lower limit cutoff for the
    890 faint object analysis would typically be S/N of 2 - 4.  In this stage,
    891 PSPhot can perform one of three types of analysis.  The difference
    892 between these options is one of speed vs detail.
    893 
    894 In the first option, PSPhot can repeat the analysis described above in
    895 sections XXX and XXX, performing a PSF fit followed by a non-PSF fit
    896 to the objects which are not PSF-like, and subtracting them.  The
    897 advantage of this option is that the faint objects are treated
    898 identically to the bright objects, and all potential parameters are
    899 measured, even for marginally extended sources.  The disadvantage of
    900 this option is that the low signal-to-noise of the objects in this
    901 stage limits the amount of information which can be extracted from
    902 them.  The marginal gain may not be worth the large expense of
    903 processing time.
    904 
    905 In the second option, PSPhot can perform only the PSF model fit to the
    906 remaining peaks, but ignore any further questions of the shape of the
    907 objects.  The advantage of this option is that it is substantially
    908 faster than performing the more complex (and less stable)
    909 multi-parameter non-linear fits on all faint objects.  On the
    910 downside, less information is learned about the objects.
    911 
    912 Finally, PSPhot can simply ignore the fitting process and instead
    913 glean information about the fainter sources on the basis of the peak
    914 characteristics.  In this option, the image is smoothed with the PSF
    915 model, and the peak for each object is measured.  The peak flux and
    916 the local peak curvature theoretically give sufficient information to
    917 recover the object flux, the centroid coordinates, and the centroid
    918 errors.  The advantage of the stage is speed, especially for the very
    919 faintest levels: if the lower limit is not sufficiently faint, the
    920 time spent in performing the smoothing (3 FFTs) cannot make up for the
    921 time which would have been spent applying the PSF model to the peaks.
    922 The downside of this method is an increased sensitivity to the local
    923 sky model (the local sky must be correctly subtracted) and fewer
    924 constraints on the quality of the detection (no Chi-Square is
    925 measured, for example).
    926 
    927 \note{In the ideal case, if we were only interested in detecting PSFs,
    928 and we had a good model for the PSF, we could optimally find the
    929 sources by smoothing the image and the variance image with the PSF model.
    930 \em write out the description of Nick's optimal PSF finding}.
    931 
    932 PSPhot allows the user to select between these three options for the
    933 analysis of the faint sources.  Three separate user-controlled
    934 signal-to-noise ratio limits are defined.  One specifies the depth to
    935 which the PSF / non-PSF analysis is performed.  A second (which must
    936 be smaller) specifies the depth to which only the PSF is fitted.  A
    937 third specifies the depth to which the analysis is performed using on
    938 the peak statistics.  If two of these are identical, then certain of
    939 these options are simply skipped.  For example, if the peak analysis
    940 threshold is set to the same value as the PSF-only threshold, no peak
    941 analysis is performed.
     1058object analysis would S/N of 5 - 10.  \note{PV3 value is 20.0?}  The
     1059lower limit cutoff for the faint object analysis would typically be
     1060S/N of 2 - 4.  \note{PV3 value is 5.0?}  Objects detected in the faint
     1061object stage are fitted with the PSF model using the linear, ensemble
     1062fitting process.
    9421063
    9431064\subsection{Aperture Correction Measurement}
    … …  
    9901111saturation. 
    9911112
    992 \note{this discussion sucks: put in some more details of my point:
    993   amplitude of systematic vs random sky errors}
    994 
    995 How important is this effect?  Consider a typical bright object with a
    996 flux of (say) 40,000 counts in an image of background 1000 counts per
    997 pixel, with FWHM of 4 pixels.  In principle, the flux of this object
    998 should be measurable with an accuracy of roughly 0.57\%
    999 ($\frac{\sqrt{40000 + 1000 \times 12}}{40000}$).  However, the
    1000 measurement of the sky is limited at some finite level by Poisson
    1001 statistics.  If we are required to use an aperture of (say) 25 pixels
    1002 in radius (eg, 5 arcseconds for an 0.2 arcsec / pixel detector), and
    1003 we have an annulus of twice this radius to measure the local sky, then
    1004 we will have an error of XXX.
    1005 
    1006 \note{outline the variation of {\em ApResid} as a function of
    1007 magnitude}.
     1113% How important is this effect?  Consider a typical bright object with a
     1114% flux of (say) 40,000 counts in an image of background 1000 counts per
     1115% pixel, with FWHM of 4 pixels.  In principle, the flux of this object
     1116% should be measurable with an accuracy of roughly 0.57\%
     1117% ($\frac{\sqrt{40000 + 1000 \times 12}}{40000}$).  However, the
     1118% measurement of the sky is limited at some finite level by Poisson
     1119% statistics.  If we are required to use an aperture of (say) 25 pixels
     1120% in radius (eg, 5 arcseconds for an 0.2 arcsec / pixel detector), and
     1121% we have an annulus of twice this radius to measure the local sky, then
     1122% we will have an error of XXX.
     1123%
     1124% \note{outline the variation of {\em ApResid} as a function of
     1125% magnitude}.
    10081126
    10091127PSPhot measures the aperture correction ({\em ApResid}) for every PSF
    … …  
    13251443using an FFT-based convolution \note{(examples?)}
    13261444
    1327 Recipe parameters which affect the PSF-convolved galaxy model fitting
    1328 process:
    1329 \begin{verbatim}
    1330 EXT_FIT_NSIGMA_CONV [9] : number of sigma
    1331 EXT_FIT_ITER
    1332 EXT_FIT_MIN_TOL
    1333 EXT_FIT_MAX_TOL
    1334 LMM_FIT_CHISQ_CONVERGENCE
    1335 LMM_FIT_GAIN_FACTOR_MODE
    1336 \end{verbatim}
    1337 
    13381445For the Exponential and DeVaucouleur fits, all parameters are fitted
    13391446in the non-linear minimization stage.  For the Sersic model, we do not
    … …  
    13681475\subsection{Convolved Radial Aperture Photometry}
    13691476
     1477For some science goals, a well-measured color of a galaxy is more
     1478important than an accurate total magnitude.  In the case of PS1, the image
     1479quality variations for stacks of different filters presents a serious
     1480challenge for the determination of precise colors.  PSPhot determines
     1481a set of PSF-matched radial aperture flux measurements in order to
     1482minimize the impact of the stack image quality variations.
     1483
     1484In PSPhotStack, the stack analysis version of PSPhot, the 5 filter
     1485images are processed together.  After the PSF models have been fitted
     1486and a best set of galaxy models have been determined, three sets of
     1487radial apertures are measured.  In the first set, the fluxes in the
     1488radial apertures are measured using the raw stack images.  The centers
     1489of the apertures for each object across the 5 filters are fixed so
     1490that the pixels represent the equivalent portions of the same galaxy
     1491for all 5 filters.  In this analysis, the best model for each object
     1492is subtracted from the image pixels for all objects excluding the
     1493object in consideration.  The 'best model' is \note{TBD}. 
     1494
     1495In addition to the raw radial apertures, the stack images are each
     1496convolved with a circular Gaussian with $\sigma$ chosen to yield an
     1497image with a typical FWHM of 6\arcsec.  The full set of radial
     1498apertures are again measured on these convolved images.  Again, the
     1499best object models are subtracted from the image for objects not being
     1500measured.  This subtraction includes the convolution to smooth the
     1501model to the effective FWHM of the convolved image.  The entire
     1502procedure is then repeated with a target FWHM of 8\arcsec. 
     1503
     1504\note{is the first convolution done with the Alard-Lupton technique?}
     1505
    13701506\subsection{Forced Photometry : PSFs}
    13711507
    13721508\subsection{Forced Photometry : galaxies}
    13731509
    1374 \subsection{Types of Object / PSF models currently available}
    1375 
    1376 \note{the discussion of the model types needs to be extended}
    1377 
    1378 \begin{itemize}
    1379 \item GAUSS  : Pure elliptical Gaussian
    1380 \item PGAUSS : polynomial approximation to a Gaussian (PGAUSS)
    1381 \item QGAUSS : power law with variable exponential term
    1382 \item SGAUSS :
    1383 \end{itemize}
    1384 
    1385 \note{discuss the stability issues with the galaxy model(s)}
    1386 
    13871510\subsection{Output Options}
    13881511
    … …  
    13951518\subsection{Trailed Sources}
    13961519
    1397 \subsection{Fixed / Known-position Sources}
    1398 
    13991520\subsection{Difference Images}
    14001521
    14011522The variance map for a difference image must be generated from the two
    1402 images use to construct the difference.  Otherwise, the low sky level
     1523images used to construct the difference.  Otherwise, the low sky level
    14031524will automatically result in inconsistent interpretation of the variance.
    14041525
    … …  
    14541575\section{Sample Tests}
    14551576
     1577\begin{verbatim}
     1578Configuration variables affecting the peak detection process:
     1579PEAKS_SMOOTH_SIGMA [2.5]   : Gaussian sigma of smoothing kernel, in pixels.
     1580PEAKS_SMOOTH_NSIGMA [2.0]  : Gaussian smoothing kernel window size in sigmas.
     1581PEAKS_NSIGMA_LIMIT [20.0]  : Detection limit on first pass (sigmas).
     1582PEAKS_NSIGMA_LIMIT_2 [5.0] : Detection limit on faint detection pass (sigmas).
     1583\end{verbatim}
     1584
     1585Recipe parameters which affect the PSF-convolved galaxy model fitting
     1586process:
     1587\begin{verbatim}
     1588EXT_FIT_NSIGMA_CONV [9] : number of sigma
     1589EXT_FIT_ITER
     1590EXT_FIT_MIN_TOL
     1591EXT_FIT_MAX_TOL
     1592LMM_FIT_CHISQ_CONVERGENCE
     1593LMM_FIT_GAIN_FACTOR_MODE
     1594\end{verbatim}
     1595
    14561596\end{document}
Note: See TracChangeset for help on using the changeset viewer.