Index: trunk/doc/psphot/psphot.tex
===================================================================
--- trunk/doc/psphot/psphot.tex	(revision 4902)
+++ trunk/doc/psphot/psphot.tex	(revision 6005)
@@ -142,4 +142,7 @@
   discussion and add in the references.}
 
+\note{Add discussion of the lessons learned from experience with previous
+  analysis programs}
+
 The Pan-STARRS IPP team decided that none of the existing packages met
 all of their needs, particularly given the very challenging goals of
@@ -156,19 +159,88 @@
 Python).
 
-\note{Add discussion of the lessons learned from experience with previous
-  analysis programs}
+\section{PSPhot Design Goals}
+
+PSPhot has a number of important requirements that it must meet, and a
+number of design goals which we believe will help to make usable in a
+wide range of circumstances.  The critical requirements of the
+Pan-STARRS IPP which drive the requirements for PSPhot:
 
 \begin{itemize}
-\item Flexible PSF model: functional form should be
-  easily modified.  
-\item PSF variation is fundamental : PSF representation should incorporate 2-D variations.  
-\item Speed fitting with accurate parameter guesses.
-\item Make good use of moment information to speed analysis.  
-\item Careful definition of PSF validity tests.  
-\item Careful analysis of aperture corrections.  
-\item Flexible non-PSF models.  
-\item Good code abstraction to simplify modification.
+\item {\bf 10 millimagnitude photometric accuracy}.  For PSPhot, this
+  implies that the measured photometry of stellar objects must be
+  substantially better than this 10 mmag since the photometry error
+  per image is combined with an error in the flat-field calibration
+  and an error in measuring the atmospheric effects.  We have set a
+  goal for PSPhot of 3mmag photometric consistency for bright stars
+  between pairs of images obtained in photometric conditions at the
+  same pointing, ie to remove sensitivity to flat-field errors.  This
+  goal splits the difference between the three main contributors and
+  still allows some leeway.  This requirement must be met for
+  well-sampled images and images with only modest undersampling.
+
+\item {\bf 10 milliarcsecond astrometric accuracy}. Relative
+  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
+  calibration at the 10mas level.  The error in the individual
+  measurements will be folded together with the errors introduced by
+  the optical system, the effects of seeing, and by the available
+  reference catalogs.  We have set a goal for PSPhot of 5mas
+  consistency between the true source postion and the measured
+  position given reasonable PSF variations under simulations.  This
+  level must be reached for images with 250 mas pixels, implying
+  PSPhot must introduce measurement errors less than 1/50th of a
+  pixel.
+
+\item {\bf processing time of 45 seconds} This requirement depends
+  strongly on the hardware organization, the amount of time spent on
+  other analysis steps, the density of stars per image, and the depth
+  for a given type of image.  For the sources at the faint limit (eg,
+  $5\sigma$), the average density of sources is expected to be roughly
+  $3\times10^5$ per square degree, while sources at the 20 $\sigma$
+  level may have densities of $\sim 5\times10^4$ per square degree.
+  Allowing 30 seconds for the PSPhot portion of the analysis, of which
+  15 is used for careful analysis of the brighter sources, 10 seconds
+  is used for PSF modeling and other overheads, and the remaining 5
+  seconds is used for the PSF fitting of the faintest source implies
+  that the detailed modelling may take roughly 3msec per source, and
+  the basic PSF fitting may be allowed 150 usec per source.
 \end{itemize}
 
+The design goals for PSPhot are chosen to make the program flexible,
+general, and able to meet the unknown usages cases future projects may
+require:
+
+\begin{itemize}
+\item {\bf Flexible PSF model} Different image sources require
+  different ways of representing the PSF.  Ideally, both analytical
+  and pixel-based versions should be possible.
+
+\item {\bf PSF spatial variation} Most images result in some spatial
+  PSF variations at a certain level.  The PSF representation should
+  naturally incorporate 2-D variations.
+
+\item {\bf Flexible non-PSF models} PSPhot must be able to represent
+  PSF-like objects as well as non-PSF sources.  It must be easy to add
+  new object models as interesting representations of sources are
+  invented.
+
+\item {\bf Clean code base} PSPhot should incorporate a high-degree of
+  abstraction and encapsulation so that changes to the code structure
+  can be performed without pulling the code apart and starting from scratch.
+
+\item {\bf PSF validity tests} PSPhot should include the ability to
+  choose different types of PSF models for diffent situations, or to
+  provide the user with methods for assessing the different PSF models.
+
+\item {\bf Careful aperture corrections} PSPhot must carefully measure
+  and correct for the photometric and astrometric trends introduced by
+  using analytical PSF models.
+
+\item {\bf User Configurable} PSPhot should allow users to change the
+  options easily and to allow different approaches to the analysis.
+
+\end{itemize}
+
 \section{PSPhot Analysis Process}
 
@@ -178,29 +250,38 @@
 
 \begin{itemize}
-\item image preparation - load data, characterize the image
+\item {\bf Image preparation} Load data, characterize the image
   background, load or construct noise and mask images.
 
-\item initial object detection - smooth, find peaks, measure basic
-  properties 
-
-\item PSF determination - select PSF candidates, perform model fits,
-  build PSF model from fits, select best PSF model class, determine
-  image aperture corrections.  
-
-\item Fit PSF objects - fit objects with PSFs, determine PSF validity,
-  subtract PSF-like objects
-
-\item Fit non-PSF objects - fit non-PSF model, select best model
-  class, subtract model
-
-\item low S/N sources - detect low-level sources, measure properties
-
-\item output - write out objects in selected format, write out
+\item {\bf Initial object detection} Smooth, find peaks, measure basic
+  properties
+
+\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),
+  select best model class, subtract model.
+
+\item {\bf Low S/N sources} Detect low-level sources, measure
+  properties (aperture or PSF)
+
+\item {\bf Aperture corrections} Measure the curve-of-growth, spatial
+  aperture variations, and background-error corrections.  
+
+\item {\bf Output} Write out objects in selected format, write out
   difference image, noise image, etc, as selected.
 \end{itemize}
 
-\note{the current configuration variables and some of the function
-  names are not very well chosen.  expect these to be modified as the
-  code base is cleaned.}
+Note that a given run of PSPhot \note{should} allow the user to
+perform any of these stages as an option.  For example, the PSF model
+may already be available from external information, in which case the
+PSF modeling stage can be skipped.  Or, when used as a library
+function, the image may have already been loaded and the mask and
+weight images constructed.  In some implementations, it may be
+possible to skip the initial object detection stage because only
+supplied sources are measured.  These are only some of the possible
+configurations.  The use of these different configurations depends on
+the source of the image, the desired detail and speed of the
+processing, and the level of accuracy desired from the analysis.
 
 \subsection{Image Preparation}
@@ -241,8 +322,7 @@
 \code{XMIN}, \code{XMAX}, \code{YMIN}, \code{YMAX}.
 
-\note{minimum valid data value is currently unimplemented}
-
-\note{mask values are currently hard-wired numbers : they should be
-given named enum values}
+\note{Mask values are currently hard-wired numbers.  We need a method
+  for user-defined mask values to be supplied.  PSLib needs to have a
+  mask registration system.}
 
 The noise image, if not supplied is constructed by default from the
@@ -267,5 +347,6 @@
 dependent on the form of this smoothing function.
 
-\note{Is this smoothing needed?}
+\note{Is this smoothing needed?  we could save time here by skipping
+it.}
 
 The local peaks in the smoothed image are found by first detecting
@@ -412,5 +493,5 @@
 represent the shape in some way}.
 
-\subsubsection{PSF Object Candidate Selection}
+\subsubsection{PSF Candidate Object Selection}
 
 The first stage of determining the PSF model for an image is to
@@ -446,5 +527,5 @@
 the image.  
 
-\subsubsection{PSF candidates object model fits}
+\subsubsection{PSF Candidate Object Model Fits}
 
 All candidate PSF objects are then fitted with the selected object
@@ -453,6 +534,6 @@
 non-linear fitting.  Non-linear fitting can be very computationally
 intensive, particularly for if the starting parameters are far from
-the minimization values.  PSPhot uses a the first and second moments
-to make a good guess for the centroid and shape parameters for the PSF
+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.  In order to minimize the impact of close neighbors, the noise
 values used in the fit are enhanced by a fraction of the deviation of
@@ -490,5 +571,292 @@
 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.
+objects in an image.  An approximate correction is measured here, with
+a more detailed correction measured after all object analysis is
+performed.  The PSF model with the best consistency of the aperture
+correction is judged to be the best model.
+
+\subsubsection{Basic Deblending}
+
+The collection of identified peaks is examined to find peaks which are
+'blended', that is, they are close enough together to make the
+analysis of one of the sources difficult if performed in isolation.
+Saturated stars also result in additional peaks which are likely to be
+invalid; it is useful to restrict a saturated star to a single primary
+position with associated neighboring peaks.
+
+The deblending process first searches for any peaks which are within
+the image cell of another peak.  All such groups are examined,
+starting with the brightest source.  An isophot is found about the
+primary peak which is at least \code{DEBLEND_SKY_NSIGMA} times the sky
+sigma above the local background and which is otherwise
+\code{DEBLEND_PEAK_FRACTION} of the primary peak central pixel flux.
+Any secondary sources which are contained within this isophot are
+considered to be blended peaks associated with the primary peak.  
+
+\subsection{Bright Source Analysis}
+
+After a PSF model has been determined, PSPhot performs the analysis of
+the bright objects in the image.  In this stage, all of the objects
+with an estimated signal to noise (based on the moments analysis)
+greater than a user-set threshold are analysed and subtracted from the
+image.  An optional successive stage then finds fainter sources and
+measures them as well (see Faint Source Analysis,
+Section~\ref{faintsources}).  In the bright source analysis stage, two
+major varients are available.  In the primary version, all objects are
+examined (in decending order of brightness) and an appropriate models
+is determined for each object which is then subtracted; in the
+alternate version, the objects are examined (in decending order of
+brightness) and the PSF-like objects subtracted first, then the
+extended objects are analysed on a second pass.
+
+\subsubsection{Fast Ensemble PSF Fitting}
+
+Before the detailed analysis of the objects 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
+PSF-like.  We also assume their position can be taken as the peak of a
+2D quadratic function fitted to the peak pixel and its surrounding 8
+pixels.  A single linear fit is used to simultaneously measure all
+source fluxes.  Since the local sky has been subtracted, this
+measurement assumes the local sky is zero.  
+
+\[
+\chi^2 = \sum_{\rm pixels} (F_{x,y} - \sum_{\rm sources} A_i PSF[x,y])^2
+\]
+
+Minimizing this equation with respect to each of the $A_i$ values
+results in a matrix equation:
+\[ M_{i,j} \bar{A_i} = \bar{F_j}\]
+where $\bar{A_i}$ is the vector of $A_i$ values, the elements of
+$M_{i,j}$ consist of the dot product of the unit-flux PSF for source
+$i$ and source $i$, and $\bar{F_j}$ is the dot product of the
+unit-flux PSF for source $i$ with the pixels corresponding to source
+$i$.  The dot products are calculated only using pixels within the
+source apertures.  Since most sources have no overlap with most other
+sources, this matrix equation results in a very sparse, mostly
+diagonal square matrix.  The dimension is the number of sources,
+likely to be 1000s or 10,000s.  Such a matrix does not lend itself to
+direct inversion.  However, an interative solution quickly yields a
+result with sufficient accuracy.  In the iterative solution, a guess
+at the solution is made; the guess is multiplied by the matrix, and
+the result compared with the observed vector $\bar{F_j}$.  The
+difference is used to modify the initial guess. This proces is
+repeated several times to achieve a good convergence.  
+
+Once a solution set for $A_i$ is found, all of the objects are
+subtracted from the by applying these values to the unit-flux PSF.
+
+\subsubsection{PSF Model applied to detected objects}
+
+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
+ratio (\note{KEYWORD}) with the PSF model.  For these fits, the
+dependent parameters are fixed by the PSF model and only the 4
+independent object model parameters are allowed to vary in the fit.
+PSPhot again uses the Levenberg-Marqardt process for the non-linear
+fitting.  The objects 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, PSPhot attempts to judge the
+quality of the PSF model as a representation of the object shape.  To
+do this, it calculates the next step of the minimization {\em allowing
+the shape parameters to vary}.  This step, essentially the
+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
+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
+elliptical Gaussian, these two parameters are $\sigma_x$ and
+$\sigma_y$.  For a generic model, the shape parameters may be defined
+differently, but the should always be two parameters which scale the
+object size in two dimensions (what about a polar-coordinate form?)
+Currently, PSPhot requires the two relevant shape parameters to be the
+first two dependent parameters in the list of model parameters (ie,
+parameters 4 \& 5).
+
+The expected distribution of the variation of the two shape parameters
+will be a function of the signal-to-noise of the object in question
+and the value of the shape parameter itself.  The expected standard
+deviation on the shape parameter is, eg, $\sigma_x / {\rm SN}$.  If
+the object 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
+step, how many standard deviations from the expected value (of 0.0) is
+this particular value?  Objects 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
+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
+as a likely defect.  
+
+At this stage of the analysis, PSPhot uses two additional indicators
+to identify good and poor PSF fits.  The first of these is the
+signal-to-noise ratio.  It is possible for the peak finding algorithm
+to identify peaks in locations which are not actually a normal peak.
+Some of these cases are in the edges of saturated, bleeding columns
+from bright stars, in the nearly flat halos of very bright stars, and
+so on.  In these cases, a local peak exists, with a lower nearby sky
+region.  However, the fitted PSF model cannot converge on the peak
+because it is very poorly defined (perhaps only existing in the
+smoothed image).  The fit can either fail to converge or it can
+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
+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
+exclusion stage are subtracted from the image.  The subtraction
+process modifies the image pixels (removing the fitted flux, though
+not the locally fitted background) but does not modify the mask or the
+noise images.  The signal-to-noise ratio in the image after
+subtraction represents the significance of the remaining flux.  If the
+subtractions are sufficiently accurate models of the PSF flux
+distribution, the remaining flux should be below 1 $\sigma$
+significance.  In practice the cores of bright stars are poorly
+represented and may have larger residual significance. \note{in future
+work, we may choose to enhance the noise to minimize detection of
+objects in the residuals of brighter objects}.
+
+\subsubsection{Blended Sources}
+
+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.
+The resulting fits are evaluated independently and any which are
+determined to be PSFs are subtracted from the image.
+
+\subsubsection{Double Sources}
+
+Sources which are judged to be non-PSF-like are confronted with two
+possible alternative choices.  First, the object 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
+otherwise not distinguished.  The initial guess for the two peaks is
+made by splitting the flux of the single source in half and locating
+the two starting peaks at +/- 2 pixels from the original peak along
+the direction of the semi-major axis of the sources, as measured from
+the second moments.  In order for the two-source model to be accepted,
+both sources must be judged as a valid PSF source.  Otherwise, the
+double-PSF model is rejected and the source is fitted with the
+available non-PSF model or models.
+
+\note{better description of the acceptance criteria; the FLT model is
+  tried before the DBL is accepted or rejected}. 
+
+\subsubsection{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
+now be fit with appropriate models for the galaxies (or other likely
+extended shapes).  Again, the fitting stage starts with the brightest
+sources (as judged by the rough S/N measured from the moments
+aperture) and working to a user defined S/N limit.  
+
+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
+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
+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
+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
+for this decision will depend on the details of the model, and so this
+level of abstraction is needed.  For example, in some case, the range
+of valid values for each of the parameters must be considered in the
+fit assessment.  In other cases, we may choose to use only the
+parameter errors and the fit Chi-Square value.
+
+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
+is subtracted from the image.  Again, the noise image is (currently)
+not modified.  
+
+\note{we have no code yet to select the best of several models for a
+  given objects.  The relative value of the Chi-Square is the obvious
+  test in this case}.
+
+\subsection{Faint Sources}
+
+\note{this is not done : should use the ensemble PSF fitting to fit
+  just the new significant peaks}
+
+After a first pass through the image, in which the brighter sources
+above a high threshold level have been detected, measured, and
+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
+noise 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.  A typical threshold for the bright
+object analysis would S/N of 5 - 10.  The lower limit cutoff for the
+faint object analysis would typically be S/N of 2 - 4.  In this stage,
+PSPhot can perform one of three types of analysis.  The difference
+between these options is one of speed vs detail.
+
+In the first option, PSPhot can repeat the analysis described above in
+sections XXX and XXX, performing a PSF fit followed by a non-PSF fit
+to the objects which are not PSF-like, and subtracting them.  The
+advantage of this option is that the faint objects are treated
+identically to the bright objects, and all potential parameters are
+measured, even for marginally extended sources.  The disadvantage of
+this option is that the low signal-to-noise of the objects in this
+stage limits the amount of information which can be extracted from
+them.  The marginal gain may not be worth the large expense of
+processing time.
+
+In the second option, PSPhot can perform only the PSF model fit to the
+remaining peaks, but ignore any further questions of the shape of the
+objects.  The advantage of this option is that it is substantially
+faster than performing the more complex (and less stable)
+multi-parameter non-linear fits on all faint objects.  On the
+downside, less information is learned about the objects.
+
+Finally, PSPhot can simply ignore the fitting process and instead
+glean information about the fainter sources on the basis of the peak
+characteristics.  In this option, the image is smoothed with the PSF
+model, and the peak for each object is measured.  The peak flux and
+the local peak curvature theoretically give sufficient information to
+recover the object flux, the centroid coordinates, and the centroid
+errors.  The advantage of the stage is speed, especially for the very
+faintest levels: if the lower limit is not sufficiently faint, the
+time spent in performing the smoothing (3 FFTs) cannot make up for the
+time which would have been spent applying the PSF model to the peaks.
+The downside of this method is an increased sensitivity to the local
+sky model (the local sky must be correctly subtracted) and fewer
+constraints on the quality of the detection (no Chi-Square is
+measured, for example).
+
+\note{In the ideal case, if we were only interested in detecting PSFs,
+and we had a good model for the PSF, we could optimally find the
+sources by smoothing the image and the noise image with the PSF model.
+\em write out the description of Nick's optimal PSF finding}.
+
+PSPhot allows the user to select between these three options for the
+analysis of the faint sources.  Three separate user-controlled
+signal-to-noise ratio limits are defined.  One specifies the depth to
+which the PSF / non-PSF analysis is performed.  A second (which must
+be smaller) specifies the depth to which only the PSF is fitted.  A
+third specifies the depth to which the analysis is performed using on
+the peak statistics.  If two of these are identical, then certain of
+these options are simply skipped.  For example, if the peak analysis
+threshold is set to the same value as the PSF-only threshold, no peak
+analysis is performed.
+
+\subsection{Aperture Correction Measurement}
 
 The important concept here is that an analytical model will always
@@ -581,80 +949,4 @@
 tested.
 
-\subsubsection{PSF Model applied to detected objects}
-
-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
-ratio (\note{KEYWORD}) with the PSF model.  For these fits, the
-dependent parameters are fixed by the PSF model and only the 4
-independent object model parameters are allowed to vary in the fit.
-PSPhot again uses the Levenberg-Marqardt process for the non-linear
-fitting.  The objects 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, PSPhot attempts to judge the
-quality of the PSF model as a representation of the object shape.  To
-do this, it calculates the next step of the minimization {\em allowing
-the shape parameters to vary}.  This step, essentially the
-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
-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
-elliptical Gaussian, these two parameters are $\sigma_x$ and
-$\sigma_y$.  For a generic model, the shape parameters may be defined
-differently, but the should always be two parameters which scale the
-object size in two dimensions (what about a polar-coordinate form?)
-Currently, PSPhot requires the two relevant shape parameters to be the
-first two dependent parameters in the list of model parameters (ie,
-parameters 4 \& 5).
-
-The expected distribution of the variation of the two shape parameters
-will be a function of the signal-to-noise of the object in question
-and the value of the shape parameter itself.  The expected standard
-deviation on the shape parameter is, eg, $\sigma_x / {\rm SN}$.  If
-the object 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
-step, how many standard deviations from the expected value (of 0.0) is
-this particular value?  Objects 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
-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
-as a likely defect.  
-
-At this stage of the analysis, PSPhot uses two additional indicators
-to identify good and poor PSF fits.  The first of these is the
-signal-to-noise ratio.  It is possible for the peak finding algorithm
-to identify peaks in locations which are not actually a normal peak.
-Some of these cases are in the edges of saturated, bleeding columns
-from bright stars, in the nearly flat halos of very bright stars, and
-so on.  In these cases, a local peak exists, with a lower nearby sky
-region.  However, the fitted PSF model cannot converge on the peak
-because it is very poorly defined (perhaps only existing in the
-smoothed image).  The fit can either fail to converge or it can
-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
-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
-exclusion stage are subtracted from the image.  The subtraction
-process modifies the image pixels (removing the fitted flux, though
-not the locally fitted background) but does not modify the mask or the
-noise images.  The signal-to-noise ratio in the image after
-subtraction represents the significance of the remaining flux.  If the
-subtractions are sufficiently accurate models of the PSF flux
-distribution, the remaining flux should be below 1 $\sigma$
-significance.  In practice the cores of bright stars are poorly
-represented and may have larger residual significance. \note{in future
-work, we may choose to enhance the noise to minimize detection of
-objects in the residuals of brighter objects}.
-
 \subsubsection{Types of Object / PSF models currently available}
 
@@ -662,118 +954,23 @@
 
 \begin{itemize}
-\item Pure elliptical Gaussian (GAUSS)
-\item polynomial approximation to a Gaussian (PGAUSS)
-\item RGAUSS
-\item QGAUSS - power law with variable exponential term
+\item GAUSS  : Pure elliptical Gaussian
+\item PGAUSS : polynomial approximation to a Gaussian (PGAUSS)
+\item QGAUSS : power law with variable exponential term
+\item SGAUSS : 
 \end{itemize}
 
-\subsection{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
-now be fit with appropriate models for the galaxies (or other likely
-extended shapes).  Again, the fitting stage starts with the brightest
-sources (as judged by the rough S/N measured from the moments
-aperture) and working to a user defined S/N limit.  
-
-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
-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
-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
-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
-for this decision will depend on the details of the model, and so this
-level of abstraction is needed.  For example, in some case, the range
-of valid values for each of the parameters must be considered in the
-fit assessment.  In other cases, we may choose to use only the
-parameter errors and the fit Chi-Square value.
-
-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
-is subtracted from the image.  Again, the noise image is (currently)
-not modified.  
-
-\note{we have no code yet to select the best of several models for a
-  given objects.  The relative value of the Chi-Square is the obvious
-  test in this case}.
-
-\subsection{Faint Sources}
-
-\note{the following discussion is theoretical : it is not yet coded}
-
-After a first pass through the image, in which the brighter sources
-above a high threshold level have been detected, measured, and
-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
-noise 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.  A typical threshold for the bright
-object analysis would S/N of 5 - 10.  The lower limit cutoff for the
-faint object analysis would typically be S/N of 2 - 4.  In this stage,
-PSPhot can perform one of three types of analysis.  The difference
-between these options is one of speed vs detail.
-
-In the first option, PSPhot can repeat the analysis described above in
-sections XXX and XXX, performing a PSF fit followed by a non-PSF fit
-to the objects which are not PSF-like, and subtracting them.  The
-advantage of this option is that the faint objects are treated
-identically to the bright objects, and all potential parameters are
-measured, even for marginally extended sources.  The disadvantage of
-this option is that the low signal-to-noise of the objects in this
-stage limits the amount of information which can be extracted from
-them.  The marginal gain may not be worth the large expense of
-processing time.
-
-In the second option, PSPhot can perform only the PSF model fit to the
-remaining peaks, but ignore any further questions of the shape of the
-objects.  The advantage of this option is that it is substantially
-faster than performing the more complex (and less stable)
-multi-parameter non-linear fits on all faint objects.  On the
-downside, less information is learned about the objects.
-
-Finally, PSPhot can simply ignore the fitting process and instead
-glean information about the fainter sources on the basis of the peak
-characteristics.  In this option, the image is smoothed with the PSF
-model, and the peak for each object is measured.  The peak flux and
-the local peak curvature theoretically give sufficient information to
-recover the object flux, the centroid coordinates, and the centroid
-errors.  The advantage of the stage is speed, especially for the very
-faintest levels: if the lower limit is not sufficiently faint, the
-time spent in performing the smoothing (3 FFTs) cannot make up for the
-time which would have been spent applying the PSF model to the peaks.
-The downside of this method is an increased sensitivity to the local
-sky model (the local sky must be correctly subtracted) and fewer
-constraints on the quality of the detection (no Chi-Square is
-measured, for example).
-
-\note{In the ideal case, if we were only interested in detecting PSFs,
-and we had a good model for the PSF, we could optimally find the
-sources by smoothing the image and the noise image with the PSF model.
-\em write out the description of Nick's optimal PSF finding}.
-
-PSPhot allows the user to select between these three options for the
-analysis of the faint sources.  Three separate user-controlled
-signal-to-noise ratio limits are defined.  One specifies the depth to
-which the PSF / non-PSF analysis is performed.  A second (which must
-be smaller) specifies the depth to which only the PSF is fitted.  A
-third specifies the depth to which the analysis is performed using on
-the peak statistics.  If two of these are identical, then certain of
-these options are simply skipped.  For example, if the peak analysis
-threshold is set to the same value as the PSF-only threshold, no peak
-analysis is performed.
+\note{discuss the stability issues with the galaxy model(s)}
+
+\subsection{Output Options}
+
+\note{need to discuss tests}
+
+\note{need to discuss failings and holes}
+
+\section{Alternative Scenarios}
+
+\subsection{Trailed Sources}
+
+\subsection{Fixed / Known-position Sources}
 
 \subsection{Difference Images}
@@ -833,7 +1030,186 @@
 discussed in the section on Galaxy models.
 
-\note{need to discuss tests}
-
-\note{need to discuss failings and holes}
+\section{PSPhot Structures and Data Elements}
+
+The following structures are described in detail in the document
+`Pan-STARRS PS-1 Image Processing Pipeline Modules Supplementary
+Design Requirements' (psModules SDRS; PSDC-430-012).
+
+\begin{datatype}
+ pmModel
+ pmModelGroup
+ pmGrowthCurve
+ pmPSF
+ pmPSFTry
+ pmSource
+ pmPeak
+ pmMoments
+\end{datatype}
+
+\note{psphot is supposed to operate on individual readouts, and use
+  the techniques used by ppImage to extract header-related metadata.
+  currently, psphot uses an alternative to the psReadout until the
+  ppImage code can be folded together with psphot}. 
+
+\subsection{Top-Level APIs}
+
+\begin{prototype}
+psMetadata     *psphotArguments (int *argc, char **argv);
+\end{prototype}
+Load the command-line arguments, parse the configuration file, and
+place the configuration information on a single metadata structure.
+This function searches for the following command line option flags,
+and places their corresponding values on the output metadata with the
+given name.  These options override any such values in the
+configuration file.
+\begin{verbatim}
+-mask (filename)      : MASK_IMAGE
+-weight (filename)    : WEIGHT_IMAGE
+-resid (filename)     : RESID_IMAGE
+-region [x0:x1,y0:y1] : ANALYSIS_REGIONP
+-photcode (code)      : PHOTCODE
+-psf (filename)       : PSF_INPUT_FILE
+-modeltest x y        : TEST_FIT_X, TEST_FIT_Y
+-model (name)         : TEST_FIT_MODEL
+-fitmode (name)       : TEST_FIT_MODE
+-fitset (name)        : TEST_FIT_SET
+\end{verbatim}
+
+The following option flags can be used to set any option:
+\begin{verbatim}
+-D  (key) (value)      : any string value
+-Df (key) (value)      : any F32 value
+-Di (key) (value)      : any S32 value
+\end{verbatim}
+
+The function next examines the remaining command-line arguments and
+complains if there are not exactly 3 arguments, reporting the program
+usage.  It sets default configuration variables, and then loads the
+configuration file specified as the third command-line option.
+Finally, it sets the \code{IMAGE} and \code{OUTPUT_FILE} config
+options to arguments 1 and 2, respecitively.
+
+\begin{prototype}
+eamReadout     *psphotSetup (psMetadata *config);
+\end{prototype}
+This function examines the configuration data in \code{config} and
+loads the image into memory.  It constructs the weight and mask images
+if they have not been specified, or loads the specified images.  The
+weight image is built based on the read noise and gain of the image,
+as extracted from the header or from the configuration options
+directly.  It defines the mask based on the selection image region,
+the values for saturation and the \code{min_VALID_PIXEL}.  
+
+\begin{prototype}
+bool            psphotModelTest (eamReadout *imdata, psMetadata *config);
+\end{prototype}
+This function is an optional test mode for psphot.  If the test mode
+has been selected, this function will attempt to fit a single object
+with the requested model.  It writes out subimage containing the
+source, the difference, the mask, and the weight.  This function may
+load a PSF model or fit a floating model.
+
+\begin{prototype}
+psStats        *psphotImageStats (eamReadout *imdata, psMetadata *config);
+\end{prototype}
+Measure the basic image properties: median sky, expected sky sigma
+
+\begin{prototype}
+psPolynomial2D *psphotImageBackground (eamReadout *imdata, psMetadata *config, psStats *sky);
+\end{prototype}
+Model the image background as a 2D polynomial and subtract from the
+image.   The should use a more sophisticated model and return the
+subtracted image.
+
+\begin{prototype}
+psArray        *psphotFindPeaks (eamReadout *imdata, psMetadata *config, psStats *sky);
+\end{prototype}
+Create a smoothed image and find all local peaks above the threshold
+level (uses: \code{PEAKS_SMOOTH_SIGMA, PEAKS_SMOOTH_NSIGMA,
+PEAKS_NSIGMA_LIMIT, PEAKS_OUTPUT_FILE})
+
+\begin{prototype}
+psArray        *psphotSourceStats (eamReadout *imdata, psMetadata *config, psArray *allpeaks);
+\end{prototype}
+Create the basic source structures for all peaks, define the initial
+pixels, measure the local sky (sky offset) and the source moments.
+
+\begin{prototype}
+bool            psphotRoughClass (psArray *sources, psMetadata *config);
+\end{prototype}
+Find the PSF clump and make the first cut source identifications
+
+\begin{prototype}
+bool            psphotBasicDeblend (psArray *sources, psMetadata *config, psStats *sky);
+\end{prototype}
+Find all blended peaks and tag, group with single primary source.
+
+\begin{prototype}
+pmPSF          *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky);
+\end{prototype}
+Try each of the selected PSF models on a subset of likely PSF stars.
+Measure the metric (aperture residual scatter) for each PSF model and
+choose the best model.  
+
+\begin{prototype}
+bool 	        psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+\end{prototype}
+Perform simultaneous fitting to all sources in the array using a
+linear fitting process which assumes all sources are PSFs and their
+positions are fixed.  Set the positions based on the bilinear
+interpolation of the peak implied by the 3x3 square of pixels
+containing the peak.  Local sky is also assumed to be correctly subtracted.
+
+\begin{prototype}
+bool            psphotFullFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+\end{prototype}
+Fit all sources in sequence starting from the brightest, and
+subtracting the sources as they are fitted.  This function only
+attempts single PSF and single EXT models and chooses between them.
+The sources are assumed to have been subtracted in advance (ie, using
+psphotEnsembleFit).  The models which do not succeed are re-subtracted
+using the prior model.
+
+\begin{prototype}
+bool            psphotBlendFit (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+\end{prototype}
+Fit all sources in sequence starting from the brightest, and
+subtracting the sources as they are fitted.  This function attempts a
+multi-source fit for blended sources, or a single PSF if it is not a
+blend, followed by both EXT and DBL models and chooses between them.
+The sources are assumed to have been subtracted in advance (ie, using
+psphotEnsembleFit).  The models which do not succeed are re-subtracted
+using the prior model.
+
+\begin{prototype}
+bool            psphotReplaceUnfit (psArray *sources);
+\end{prototype}
+After models have been attempted for all sources, this function
+replaces the sources which were temporarily subtracted, but which did
+not succeed or converge on a good solution.
+
+\begin{prototype}
+bool            psphotApplyPSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+\end{prototype}
+Attempt to fit the PSF model to all sources in brightness order,
+subtracting the resulting model if successful.  Only attempts single
+PSF models.   
+
+\begin{prototype}
+bool            psphotFitExtended (eamReadout *imdata, psMetadata *config, psArray *sources, psStats *skyStats);
+\end{prototype}
+Attempt to fit the PSF model to all sources in brightness order,
+subtracting the resulting model if successful.  Only attempts single
+EXT models.
+
+\begin{prototype}
+bool            psphotApResid (eamReadout *imdata, psArray *sources, psMetadata *config, pmPSF *psf);
+ \end{prototype}
+Measure the curve-of-growth and the aperture correction trend.
+
+\begin{prototype}
+void            psphotOutput (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+\end{prototype}
+Write out data in various formats as selected.
 
 \section{User's Guide}
@@ -879,3 +1255,13 @@
 \section{Sample Tests}
 
+\section{Further Work to be Completed}
+
+\begin{itemize}
+\item convert to pmCell as input data
+\item loop over all readouts in a pmCell
+\item write out multiple files?
+\item better method for defining the recipe?
+\item additional options for image background
+\item image background should return a background image
+
 \end{document}
