Index: /trunk/doc/misc/Objects.tex
===================================================================
--- /trunk/doc/misc/Objects.tex	(revision 2940)
+++ /trunk/doc/misc/Objects.tex	(revision 2940)
@@ -0,0 +1,277 @@
+\documentclass[panstarrs,spec]{panstarrs}
+
+\title{PSPhot} % put in your title
+\subtitle{A Strategy for Astronomical Object Detection, Classification, and Measurement}
+\author{Eugene Magnier}
+\audience{Pan-STARRS Science}
+\shorttitle{PSPhot}
+\group{Pan-STARRS IPP}
+\project{Pan-STARRS IPP}
+\organization{Institute for Astronomy}
+\version{DR}
+\docnumber{PSDC-xxx-xxx}
+
+\begin{document}
+\maketitle
+
+\section{}
+
+Object Detection
+
+Object detection consists of smoothing the image with a kernel
+approximately equivalent to the PSF and detecting the significant
+peaks above a threshold.  The detection efficiency is only
+significantly affected by the accuracy of the PSF model in this stage
+at the faint end of the detection process.  For brighter sources, a
+very rough approximation will do.  I have implemented this step by
+smoothing with a 1D gaussian first in the X direction then in the Y
+direction.  With a gaussian with $\sigma = 1.0$ pixels, and extending
+only to $2 \sigma$ from the center, this process takes roughly 1.5
+seconds on {\em alala} for a Megacam CCD, equivalent to roughly 2.8
+seconds for a full Gigapixel camera 4k x 4k chip.  
+
+My peak-finding routine runs first along each row and finds all peaks
+above a requested threshold in the row data.  It then examines each
+row peak and eliminates any peaks which are not local peaks when all 8
+neighbor pixels are considered.  If a peak has any neighbors which are
+at the same flux level, it is kept.  This last point is a weakness of
+the algorithm as it leaves behind multiple entries for saturated
+regions, and it does not attempt to exclude double peaks from
+individual stars (typical for slightly defocussed images).  The peak
+finding and cleaning step takes a total of roughly 25 msec per
+resulting output peak.  This analysis should scale like the number of
+peaks for large numbers of peaks or like the number of pixels for
+sparse images with few peaks.  
+
+Object Classification
+
+After objects have been detected by finding their peaks, some basic
+classification should be performed.  This will allow us to identify a
+subset of objects which are likely to be stars, those which are likely
+to be cosmic rays and defects, those which are likely to be galaxies,
+and those which are likely to be saturated regions.  This
+classification can be acheived on the basis of simple object moments
+measured for a small aperture centered on the object.  My
+implementation of the object moment analysis uses a circular aperture
+(of 5 pixels) centered on the peak and an external annulus (of 15 - 25
+pixels) to measure the 'sky' or background level.  Alternatively, the
+sky could be measured once for the entire image in large (eg, 64x64
+pixel) boxes.  This latter choice is faster only if the image contains
+a large number of objects.  With 1256 pixels in the sky annulus, the
+break-even point between these two options is for $\sim 7300$ objects
+for a megacam image or $\sim 13000$ objects for a full GPC image,
+given the same aperture sizes.  Note also that the GPC images from
+Pan-STARRS are likely to have much smaller PSF sizes compared with
+Megacam (in terms of pixels: typical FWHM is expected to be $\sim$ 2.5
+vs 3.8 pixel for Megacam).  
+
+Objects can be well-classified in an initial pass simply on the basis
+of their moments and other simple statistics.  In this pass, the goal
+is to identify a collection of likely stars which can be used to
+measure the image PSF model.  This requirement suggests several
+natural filters: we should select objects for which the total number
+of counts and sky signal imply a signal-to-noise ratio above some
+generous threshold.  For my analysis, I have used $S/N > 50$.  We
+should also reject objects with saturated pixels: those with peak
+pixel values above the saturation limit.  I rejected all objects with
+peak flux $> 60000$.  We should also avoid objects which are too close
+to the edges.  Most significantly, we should select objects which have
+second moments which suggest they are stellar.  Since we assume that
+the stellar objects have roughly consistent shape (ie, PSF variations
+are less than some limit), then the stellar objects should all fall
+within a narrow region in the second-moment space.  In my analysis, I
+construct a grid of $\sigma_x, \sigma_y$ (second moments in the x and
+y directions) and count the number of objects in each grid box.  The
+peak box is assumed to represent the stellar objects.  I then measure
+the scatter of both $\sigma_x, \sigma_y$ for all objects within a
+certain range of the peak value (0.2, 0.2 pixels).  I then select only
+objects which like within $\pm 2\sigma$ of the median of this
+collection.  The resulting collection of objects is a very likely
+sample of stellar objects, which are good candidates from which to
+measure the image PSF model.
+
+We can also roughly classify the other objects based on this
+second-moment distribution space.  Objects with $\sigma_x <<
+\bar{\sigma_x}$ or $\sigma_y << \bar{\sigma_y}$ are likely to be
+cosmic-rays or other instrumental defects.  Objects with $\sigma_x >>
+\bar{\sigma_x}$ or $\sigma_y >> \bar{\sigma_y}$ are likely to be
+galaxies.  Other object types are difficult to distinguish from the
+simple moments, and must be tested on a case-by-case basis.
+
+Object Measurements
+
+We have several goals in measuring the objects in the image.  We would
+like to improve our determination of the object type.  Several basic
+object types may be identified: 1) objects which are well-represented
+with a stellar PSF, 2) objects which are sub-PSF in some dimension and
+thus must be a defect, 3) objects which are well-represented with one
+or more galaxy flux models, 4) objects which represent other
+recognizable artifacts such as highly saturated stars, bleeding
+pcolumns, satelite trails, etc.
+
+For stellar objects, we need to measure the object photometry and its
+astrometry.  An important trade-off to be made here is between
+aperture photometry and PSF modeling.  Aperture photometry has the
+advantage of exactly adding up the light from the star of interest. In
+principal, if the field is empty of other stars, and a large enough
+aperture is used, and the signal-to-noise in the star is very high,
+aperture photometry will be limited only by the photon statistics of
+the star in question.  In practice, several factors of varying
+importance break this ideal: There are always neighboring objects
+which contaminate the measurement.  For faint objects, the aperture
+must be small to maximize the signal-to-noise contribution relative to
+the sky, entailing aperture corrections.  PSF modeling is less
+sensitive to the presence of neighbors (they may be naturally
+excluded) or the choice of sky aperture.  It is also more robust for
+the photometry of faint objects.  It has several disadvantages as
+well: the PSF model is never a perfect representation of the PSF,
+entailing aperture corrections.  The more accurately the PSF model
+represents the actual stellar profile, in general, the more
+time-consuming will be the fitting process.  
+
+I have investigated several representations of the PSF as an
+analytical function.  A traditional technique is to fit the objects
+with a 2-D Gaussian function:
+\[
+PSF(x,y) = Z_{\rm pk} e^{-z} + Sky
+\]
+where
+\[
+z = \frac{x^2}{2\sigma_x^2} + \frac{y^2}{2\sigma_y^2} + x y \sigma_{xy}
+\]
+where $x = X - X_o$ and $y = y - Y_o$.  The old-standby program,
+DoPhot, uses, as a default representation, a polynomial approximation
+to the Gaussian by defining:
+\[
+PSF(x,y) = Z_{\rm pk} (1 + z + z^2/2 + z^3/6)^{-1} + Sky
+\]
+with all other definitions the same as above.  This function may be
+modified by adding coefficients to the second and third order $z$
+terms to produce a function which is more or less 'wingy' as needed:
+\[
+PSF(x,y) = Z_{\rm pk} (1 + z + \beta_1 z^2/2 + \beta_2 z^3/6)^{-1} + Sky
+\]
+Several points may be noted.  First, the shift from a pure Gaussian to
+the polynomial expansion is motivated by processing speed as well as
+by the added ability to modify the shape with $\beta_1$ and $\beta_2$.
+Care must be taken, however, if the full generalization is used: there
+is a strong degeneracy between $\beta_1$ and $\beta_2$ which can drive
+the fitting routine into bad regimes.  It is usually safer to re-write
+the function as:
+\[
+PSF(x,y) = Z_{\rm pk} (1 + z + \beta_1 z^2/2 (1 + \beta_2 z/3))^{-1} + Sky
+\]
+in which case the fitting process is stabilized.  
+
+In practice, these Gaussian and pseudo-Gaussian functions do a poor
+job of representing real astronomical objects.  One useful measure of
+how well the objects are fitted by the function is the difference
+between the photometry of the fitted function and the equivalent
+aperture.  Both the amplitude of the residual (expressed as a
+magnitude) and the scatter of the residual for an image are useful
+indicators of the accuracy of the fit, and therefore the intrinsic
+photometry.  For each image, it is necessary to measure the amplitude
+of the residual in order to correct the fit photometry to a consistent
+absolute photometry system which is consistent from image to image.
+The error on this aperture correction, as represented by the scatter,
+provides a limit to the photometry which may be achieved from the
+fitted function.  The pseudo-Gaussians typically have aperture
+corrections in range of 3-6\%, with internal scatter in the vicinity
+of 1-1.5\%.  Careful selection of the stars used to measure the
+aperture residuals may reduce the amplitude of the scatter.  The
+reason for the large residuals is quite clear when one examines actual
+stellar profiles.  The stellar profiles are typically a single power
+law slope modified by a central core.  The Gaussian has an exponential
+fall-off while the pseudo-Gaussians have multiple combined power
+laws.  
+
+In addition to only poorly representing the PSF, these simple
+pseudo-Gaussians fall down because real stellar profiles frequently
+have additional structure.  In particular, it is typical that the
+inner core of the light profile, which is dominated by the effects of
+seeing is elongated in a different way (different ellipticity,
+different orientation) from the outer wings, which are dominated by
+the effects of scattering from dust particles on the mirror or in the
+atmosphere.  In order to represent these two types of structures, it
+is necessary to have two independent radial profiles.  I have
+experimented with several representations, and have found the
+following functional form to represent real images well and to be
+generally stable:
+\[
+PSF(x,y) = Z_{\rm pk} (1 + z_1^N + z_2^M)^{-1} + Sky
+\]
+where
+\[
+z_1 = \frac{x^2}{2\sigma_{x1}^2} + \frac{y^2}{2\sigma_{y1}^2} + x y \sigma_{xy1}
+z_2 = \frac{x^2}{2\sigma_{x2}^2} + \frac{y^2}{2\sigma_{y2}^2} + x y \sigma_{xy2}
+\]
+and the values of $N$ and $M$ depend on the details of the images.
+These could potentially be fitted as well with the rest of the terms,
+but such a fitting process is fairly unstable.  In a collection of
+Megacam images, I have found the values of $N = 1$ and $M = 2.25$ to
+work quite well.  Perhaps the value of $M$ could be determined for
+each image by fitting the 1-D function of $\log (f) vs \log(r)$.  
+
+Experiments with 40 Megacam images shows that the value of the
+aperture residual is much smaller and more consistent with the
+2-component model above: GIVE dAp, Sigma dAp values.
+
+Problems with the above function: UNSTABLE for some images:
+
+Poor images (bad guiding, bad image quality -> unstable fit!)
+
+Aperture Photometry vs PSF Photometry
+
+PSF Models
+
+Speed vs Detail
+
+
+PSF Model    Typical Execution Time (msec / object)
+Pgauss       1.5
+pgauss       2.3
+Sgauss       3.2
+sgauss       5.4
+Moments      0.2
+
+These numbers are based on 900 image pixels for the models, along with
+a 5 pixel object radius and a 15-25 pixel sky annulus for moments.
+The size of the analysis box should be large enough that the stellar
+contribution falls well below the sky noise level at the outer
+regions.  For some of the brighter objects (peak $>$ 10000 counts),
+this choice may be somewhat small.  For objects with peak less than
+1000 counts, the analysis box can probably be somewhat reduced.  In
+addition, the actual number of pixels used in the analysis can be
+reasonably reduced by a significant factor (2-4) by binning together
+pixels at large radii.  The slope of the stellar flux profile at large
+radii is small, so binning together pixels can increase the
+signal-to-noise per-pixel and speed up the analysis without
+diminishing the accuracy of the measurement.
+
+Stellar Densities
+
+APIs
+
+typedef struct {
+  int x;
+  int y;
+  int class;  /* an enum? */
+  float counts;
+} Peaks;
+
+class: isolated, edge, flattop
+
+typedef struct {
+  float x;
+  float y;
+  float sx;
+  float sy;
+  float sxy;
+  float zsum;
+  float zpeak;
+  float sky;
+} ObjectStats;
+
+pmFindImagePeaks (psImage image, 
+
+\end{document}