\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.  

\section{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.

\section{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
significantly diminishing the accuracy of the measurement.

Stellar Densities

Galaxy Models

Names:

Peak, Source, Object, Star, Galaxy, PSF 

Peak: a single pixel higher than its 8 neighbors

Source: 

Models:
 real Gaussian:   f(z) = Zo exp (-z) + So
 pseudo Gaussian: f(z) = Zo (1 + z + (1/2) z^2 + (1/6) z^3)^(-1)
 Waussian (?):    f(z) = Zo (1 + z + B_1 (1/2) z^2 (1 + B_2 (1/3) z))^(-1)
 Core Powerlaw:   f(z) = Zo (1 + z_1 + z_2^n)^(-1)
 moffat func:
 sersic           f(z) = Zo exp (-z^n) (n integer? n/2 integer?)
 sersic w/ core   f(z) = Zb (1 + z + (1/2) z^2 + (1/6) z^3)^(-1) + Zd exp (-A z^n)
 (n = 1/2 : exponential disk; n = 1/8 : deVau)

\end{document}
