Index: /trunk/doc/modules/ModulesSDRS.tex
===================================================================
--- /trunk/doc/modules/ModulesSDRS.tex	(revision 3067)
+++ /trunk/doc/modules/ModulesSDRS.tex	(revision 3068)
@@ -1,3 +1,3 @@
-%%% $Id: ModulesSDRS.tex,v 1.26 2004-12-24 02:30:50 price Exp $
+%%% $Id: ModulesSDRS.tex,v 1.27 2005-01-20 07:46:41 eugene Exp $
 \documentclass[panstarrs]{panstarrs}
 
@@ -1036,4 +1036,323 @@
 of the same type.
 
+\section{Object Detection, Measurement, and Classification Routines}
+
+\subsection{Overview}
+
+The process of finding, measurering, and classifying astronomical
+sources on images is one of the critical tasks of the IPP or any
+astronomical software system.  In this section, we define structures
+and functions related to the task of source detection and measurement.
+The elements defined in this section are generally low-level
+components which can be connected together to construct a complete
+object measurement suite.  An example pseudo-C program using these
+functions is provided in Appendix~\ref{psphot}.  
+
+We first define the collection of structures needed to carry
+information about the detected sources.  A major challenge is to
+define what we mean by an astronomical object in the context of image
+source detection.  An astronomical object may be as simple as a
+stellar point source, or it may consist of a more complex galaxy which
+has smooth extended structure, or it may consist of an irregular
+galaxy or galaxy group with substantial and complex sub-structure, or
+it may consist of complex non-stellar structures such as planetary
+nebular, reflection nebular, outflows and jets.  
+
+The simplest objects (ie, stars) can be sufficiently modelled by the
+point-source function (PSF).  More complex objects (such as simple,
+smooth galaxies), may have approximate analytical models which
+represent their morphology with more-or-less accuracy.  In the extreme
+cases, the objects are not well modelled at all and must be
+represented in other ways.  Thus, one aspect of our data structures
+must be elements to specify if an object has been represented by a
+model, what the model parameters are, and how well it is represented
+by the model.  Another aspect of the data structures must be a
+representation of the pixels associated with the object so complex
+structures may be reference without attempting to supply an analytical
+model.  Finally, it is often useful to allow a single complex model to
+be represented as a collection of simpler contained structures which
+may be modelled.  Thus, the representation of an object must be
+capable of identifying children, or substructures, of that object.
+
+Two additional aspects must be considered.  First, source detection
+need no be performed on a single image in isolation: it is necessary
+for multiple realizations of the same source in multiple images to be
+measured together (whether or not through simultaneous fitting in
+multiple bands or via application of the results from one image to
+another image).  Second, it may be necessary for object measurements
+to be performed on pixels in which no source is actually detectable.
+For example, this is a convenient way to provide flux upper limits at
+the locations of known objects.
+
+\subsection{Source Structures}
+
+We start by defining a single source detected in a single band:
+\begin{verbatim}
+typedef struct {
+  psPeak *peak;
+  psPixels *pixels;
+  psMoments *moments;
+  psModel *models;
+  psSourceType type;
+} psSource;
+\end{verbatim}
+
+This source has the capacity for several types of measurements.  The
+simplest measurement of a source is the location of the peak pixel
+associated with the source:
+\begin{verbatim}
+typedef struct {
+  int x;
+  int y;
+  float counts;
+  psPeakType class;
+} psPeak;
+\end{verbatim}
+
+A peak pixel may have several features which may be determined when
+the peak is found or measured.  These are carried by the
+\code{psPeakType} structure.  The \code{PS_PEAK_LONE} represents a
+single pixel which is higher than its 8 immediate neighbors.  The
+\code{PS_PEAK_EDGE} represents a peak pixel which touching the image
+edge.  The \code{PS_PEAK_FLAT} represents a peak pixel which has more
+than a specific number of neighbors at the same value, within some
+tolarence:
+\begin{verbatim}
+typedef enum {
+  PS_PEAK_LONE;
+  PS_PEAK_EDGE;
+  PS_PEAK_FLAT;
+} psPeakType; 
+\end{verbatim}
+
+In addition, the pixels which contain the source may be specified with
+the \code{psPixels} element.  A \code{psPixels} collection identifies
+a group of pixels by specifying sets of pixel spans:
+\begin{verbatim}
+typedef struct {
+  int nspan;
+  int rmin, rmax;
+  int cmin, cmax;
+  psSpan *s;
+} psPixels;
+\end{verbatim}
+%
+\begin{verbatim}
+typedef struct {
+  short row;
+  short col0;
+  short col1;
+} psSpan;
+\end{verbatim}
+
+One of the simplest measurements which can be made quickly for an
+object are the object moments.  We specify a structure to carry the
+moment information for a specific source: 
+
+\begin{verbatim}
+typedef struct {
+  float x;
+  float y;
+  float Sx;
+  float Sy;
+  float Sxy;
+  float Sum;
+  float Peak;
+  float Sky;
+  int   nPixels;
+} psMoments;
+\end{verbatim}
+
+An object may be modelled with some analytical function.  The result
+of the model includes the model parameters and their errors, along
+with the fit $\Chi^2$.
+
+\begin{verbatim}
+typedef struct {
+  psObjectModel type;
+  int Nparams;
+  float *params;
+  float *dparams;
+  float chisq;
+} psModel;
+\end{verbatim}
+
+\begin{verbatim}
+typedef enum {
+  PS_MODEL_GAUSS;
+  PS_MODEL_PGAUSS;
+  PS_MODEL_TWIST_GAUSS;
+  PS_MODEL_WAUSS;
+  PS_MODEL_SERSIC;
+  PS_MODEL_SERSIC_BULGE;
+} psModelType; 
+\end{verbatim}
+
+\subsection{Basic Object Detection APIs}
+
+\begin{verbatim}
+psVector *pmFindVectorPeaks (psVector vector, float threshold);
+\end{verbatim}
+
+Find all local peaks in the given vector above the given threshold.  A
+peak is defined as any element with a value greater than its two
+neighbors with a value above the threshold.  Two types of special
+cases have to be addressed.  If an element has the same value as the
+following element, it is not considered a peak.  If an element has the
+same value as the preceeding element (but not the following), then it
+is considered a peak.  Note that this rule (arbitrarily) identifies
+flat regions by their trailing edge.  At the edges, the element must
+be higher than its neighbor, or equal for the end edge.  This rule
+again places the peak associated with a flat region which touches the
+image edge at the image edge.  The result fo this function is a vector
+containing the coordinates of the detected peaks. 
+
+\begin{verbatim}
+psList *pmFindImagePeaks (psImage image, float threshold);
+\end{verbatim}
+
+Find all local peaks in the given image above the given threshold.
+This function should find all row peaks using
+\code{pmFindVectorPeaks}, then test each row peak and exclude peaks
+which are not local peaks.  A peak is a local peak if it is higher
+than all 8 neighbors.  \tbd{flat region?}  The result of this function
+is a list of \code{psPeak} entries.
+
+\begin{verbatim}
+psList *pmCullPeaks (psList peaks, float maxvalue, psRegion valid);
+\end{verbatim}
+
+Remove certain types of peaks from a list of peaks based on the given
+criteria.  Peaks should be eliminated if they have a peak value above
+the given maximum value limit or if the fall outside the valid
+region.  
+
+\begin{verbatim}
+psSource *pmSourceLocalSky (psImage *image, psPeak *peak, float inner_radius, float outer_radius), 
+\end{verbatim}
+
+Measure the local sky in the vicinity of the given \code{peak}.  The
+image pixels in the square annulus with inner and outer radii as
+specified are used to measure the median (clipped mean?)  (statistic
+from \code{psStats}?) in the vicinity of the the specified peak
+coordinates.  The resulting sky is applied to the \code{psMoments}
+element of the allocated \code{psSource} structure.  The input
+\code{peak} is also added to the \code{psSource} element, which is
+returned.
+
+\begin{verbatim}
+psSource *pmSourceMoments (psImage *image, psSource *source, float radius), 
+\end{verbatim}
+
+Measure source moments for the given \code{source}, using the value of
+\code{source.moments.sky} provided.  The resulting moment values are
+applied to the \code{source.moments} entry, and the source is
+retained.  The moments are measured within the given radius of the
+\code{source.peak} coordinates.
+
+\begin{verbatim}
+psSource *pmSourceRoughClass (psSource *source, float saturate, float SNlim, psRegion valid);
+\end{verbatim}
+
+\begin{verbatim}
+pmSourceFitModel (psImage *image, psSource *source, psModelType model);
+\end{verbatim}
+
+\begin{verbatim}
+pmSourceFitModel (psImage *image, psSource *source, psModelType model);
+\end{verbatim}
+
+\subsection{Basic Object Models}
+
+float PGaussian (psVector *deriv, psVector *params, psVector *x);
+
+  X = x[0] - param[2];
+  Y = x[1] - param[3];
+  
+  px = param[4]*X;
+  py = param[5]*Y;
+
+  z = 0.5*SQ(px) + 0.5*SQ(py) + param[6]*X*Y;
+  r = 1.0 / (1 + z + 0.5*z*z*(1 + z/3)); /* ~ exp (-Z) */
+  f = param[1]*r + param[0];
+  q = param[1]*r*r*(1 + z + 0.5*z*z);
+  /* note difference from gaussian: q = param[1]*r */
+
+  deriv[0] = +1;
+  deriv[1] = +r;
+  deriv[2] = q*(2*px*param[4] + param[6]*Y);
+  deriv[3] = q*(2*py*param[5] + param[6]*X);
+  deriv[4] = -2*q*px*X;
+  deriv[5] = -2*q*py*Y;
+  deriv[6] = -q*X*Y;
+
+float Waussian (psVector *deriv, psVector *params, psVector *x);
+
+  X = x - param[2];
+  Y = y - param[2];
+  
+  px = param[4]*X;
+  py = param[5]*Y;
+
+  z = 0.5*SQ(px) + 0.5*SQ(py) + param[6]*X*Y;
+  k = 0.5*z*z*(1 + param[8]*z/3);
+  r = 1.0 / (1 + z + param[7]*k); /* ~ exp (-Z) */
+  f = param[1]*r + param[0];
+  q = param[1]*r*r*(1 + param[7]*z*(1 + param[8]*z/2));
+  /* note difference from gaussian: q = param[1]*r */
+
+  deriv[0] = +1;
+  deriv[1] = +r;
+  deriv[2] = q*(2*px*param[4] + param[6]*Y);
+  deriv[3] = q*(2*py*param[5] + param[6]*X);
+  deriv[4] = -2*q*px*X;
+  deriv[5] = -2*q*py*Y;
+  deriv[6] = -q*X*Y;
+  deriv[7] = -100*param[1]*r*r*k;
+  deriv[8] = -100*param[1]*r*r*param[7]*(z*z*z)/6;
+  /* the values of 100 dampen the swing of param[7,8] */
+
+float TwistGaussian ()
+
+  X = x - param[2];
+  Y = y - param[3];
+  
+  px1 = param[4]*X;
+  py1 = param[5]*Y;
+  px2 = param[7]*X;
+  py2 = param[8]*Y;
+
+  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + param[4]*X*Y;
+  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + param[9]*X*Y;
+
+  r = 1.0 / (1 + z1 + pow(z2,Npow));
+  f = param[5]*r + param[6];
+
+  q1 = param[5]*SQ(r);
+  q2 = param[5]*SQ(r)*Npow*pow(z2,(Npow-1));
+
+  deriv[0] = +1;
+  deriv[1] = +r;
+  deriv[2] = q1*(2*px1*param[4] + param[6]*Y) + q2*(2*px2*param[7] + param[9]*Y);
+  deriv[3] = q1*(2*py1*param[5] + param[6]*X) + q2*(2*py2*param[8] + param[9]*X);
+
+  /* these fudge factors impede the growth of param[4] beyond param[7] */
+  f1 = fabs(param[7]) / fabs(param[4]);
+  f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1;
+  deriv[4] = -2*q1*px1*X*f2;
+
+  /* these fudge factors impede the growth of param[5] beyond param[8] */
+  f1 = fabs(param[8]) / fabs(param[5]);
+  f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1;
+  deriv[5] = -2*q1*py1*Y*f2;
+
+  deriv[6] = -q1*X*Y;
+
+  deriv[7] = -2*q2*px2*X;
+  deriv[8] = -2*q2*py2*Y;
+  deriv[9] = -q2*X*Y;
+
+
+
 \section{Revision Change Log}
 \input{ChangeLogSDRS.tex}