Changeset 3070 for trunk/doc/modules/ModulesSDRS.tex

Timestamp:

Jan 21, 2005, 3:59:10 PM (21 years ago)

Author:

eugene

Message:

updates for cycle 5

File:

• trunk/doc/modules/ModulesSDRS.tex (modified) (13 diffs)

trunk/doc/modules/ModulesSDRS.tex

@@ -1 @@
-%%% $Id: ModulesSDRS.tex,v 1.27 2005-01-20 07:46:41 eugene Exp $
+%%% $Id: ModulesSDRS.tex,v 1.28 2005-01-22 01:59:10 eugene Exp $
 \documentclass[panstarrs]{panstarrs}
 
@@ -11 +11 @@
 \project{Pan-STARRS Image Processing Pipeline}
 \organization{Institute for Astronomy}
-\version{02}
+\version{03}
 \docnumber{PSDC-430-012}
 
@@ -27 +27 @@
 % version  Date            Description
 DR & 2004 Jun 7 & Draft \\ \hline
-00 & 2004 Aug 16 & release for cycle 3 \\ \hline
+00 & 2004 Aug 16 & final for cycle 3 \\ \hline
 01 & 2004 Oct 12 & draft for cycle 4 \\ \hline
+02 & 2004 Nov 30 & final for cycle 4 \\ \hline
+03 & 2005 Jan 21 & draft for cycle 5 \\ \hline
 \RevisionsEnd
 
@@ -1040 +1042 @@
 \subsection{Overview}
 
-The process of finding, measurering, and classifying astronomical
+The process of finding, measuring, 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
@@ -1047 +1049 @@
 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}.  
+functions is provided in Appendix~\ref{psphot}.
 
 We first define the collection of structures needed to carry
@@ -1053 +1055 @@
 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
+stellar point source, or it may consist of a galaxy which has smooth
+extended structure; 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 nebulae, reflection
+nebulae, outflows and jets.
+
+The simplest objects (ie, stars) can be sufficiently modeled 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.
+cases, the objects are not well modeled 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
+referenced 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
+need not 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}
+another image).  Second, it will be necessary to performed object
+measurements on pixels in which no source is actually detected.  For
+example, this is a convenient way to provide flux upper limits at the
+locations of known objects.
+
+\subsection{Structures to Describe Sources}
 
 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;
+  psPeak *peak;            // description of peak pixel
+  psImage *pixels;         // rectangular region including object pixels
+  psImage *mask;           // mask to mark pixels associated with object in region
+  psMoments *moments;      // basic moments measure for the object
+  psModel *models;         // model parameters and type
+  psSourceType type;       // best identification of object
 } 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:
+simplest measurement of a source is the location and flux of the peak
+pixel associated with the source:
 \begin{verbatim}
 typedef struct {
-  int x;
-  int y;
-  float counts;
-  psPeakType class;
+  int x;                   // x-coordinate of peak pixel
+  int y;                   // y-coordinate of peak pixel
+  float counts;            // value of peak pixel (above sky?)
+  psPeakType class;        // description of peak
 } 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
+the peak is found or measured.  These are specified by the
+\code{psPeakType} enum.  \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
@@ -1120 +1123 @@
 \begin{verbatim}
 typedef enum {
-  PS_PEAK_LONE;
-  PS_PEAK_EDGE;
-  PS_PEAK_FLAT;
+  PS_PEAK_LONE;             // isolated peak
+  PS_PEAK_EDGE;             // peak on edge
+  PS_PEAK_FLAT;             // peak has equal-value neighbors
 } 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}
+The pixels which contain the source may be specified with the
+\code{psImage *pixels} element, and the mask image may be used to
+exclude any pixels which are not considered part of the source.  Note
+that the source image may be simply a subimage of the main image or a
+separate copy if the pixels are modified (eg, by subtracting flux from
+other sources, etc).
 
 One of the simplest measurements which can be made quickly for an
@@ -1152 +1142 @@
 \begin{verbatim}
 typedef struct {
-  float x;
-  float y;
-  float Sx;
-  float Sy;
-  float Sxy;
-  float Sum;
-  float Peak;
-  float Sky;
-  int   nPixels;
+  float x;                  // x-coord of centroid
+  float y;                  // y-coord of centroid
+  float Sx;                 // x-second moment
+  float Sy;                 // y-second moment 
+  float Sxy;                // xy cross moment
+  float Sum;                // pixel sum above sky (background)
+  float Peak;               // peak counts above sky
+  float Sky;                // sky level (background)
+  int   nPixels;            // number of pixels used
 } 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$.
+An object's flux distribution may be modelled with some analytical
+function.  The description of the model includes the model parameters
+and their errors, along with the fit $\chi^2$.  The model type is
+specified by the enum \code{psObjectModel}, specified below.  We
+discuss the details of these models in section~\ref{ObjectModels}.
 
 \begin{verbatim}
 typedef struct {
-  psObjectModel type;
-  int Nparams;
-  float *params;
-  float *dparams;
-  float chisq;
+  psObjectModel type;       // model to be used
+  int Nparams;              // number of parameters
+  float *params;            // parameter values
+  float *dparams;           // parameter errors
+  float chisq;              // fit chisq
 } psModel;
 \end{verbatim}
@@ -1185 +1177 @@
   PS_MODEL_WAUSS;
   PS_MODEL_SERSIC;
-  PS_MODEL_SERSIC_BULGE;
+  PS_MODEL_SERSIC_CORE;
 } psModelType; 
 \end{verbatim}
@@ -1191 +1183 @@
 \subsection{Basic Object Detection APIs}
 
-\begin{verbatim}
-psVector *pmFindVectorPeaks (psVector vector, float threshold);
+In this section, we specify a collection of basic functions which
+operate on images and sources.  We define them roughly in order in
+which we expect to use them in a basic object detection process.
+
+\begin{verbatim}
+psVector *pmFindVectorPeaks(const 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);
+neighbors and with a value above the threshold.  Two types of special
+cases must be addressed.  Equal value elements: 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.  Edge
+cases: At start of the vector, the element must be higher than its
+neighbor.  At the end of the vector, the element must be higher or
+equal to its neighbor.  These two rules again places the peak
+associated with a flat region which touches the image edge at the
+image edge.  The result of this function is a vector containing the
+coordinates (element number) of the detected peaks (type
+\code{psU32}).
+
+\begin{verbatim}
+psList *pmFindImagePeaks(const psImage *image, float threshold);
 \end{verbatim}
 
@@ -1215 +1214 @@
 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);
+which are not local peaks.  A peak is a local peak if it has a higher value
+than all 8 neighbors.  If the peak has the same value as its +y
+neighbor or +x neighbor, it is NOT a local peak.  If any other
+neighbors have an equal value, the peak is considered a valid peak.
+Note two points: first, the +x neighbor condition is already enforced
+by \code{pmFindVectorPeaks}.  Second, these rules have the effect of
+making flat-topped regions have single peaks at the (+x,+y) corner.
+When selecting the peaks, their type must also be set.  The result of
+this function is a list of \code{psPeak} entries.
+
+\tbd{do we need a function psVector *psImageRowVector (psImage *image, int row);}
+
+\begin{verbatim}
+psList *pmCullPeaks(psList *peaks, float maxvalue, const psRegion *valid);
 \end{verbatim}
 
@@ -1226 +1233 @@
 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), 
+region.  The result of the function is to reduce the number of peaks
+in the input peaks list, and return the resulting list. 
+
+\begin{verbatim}
+psSource *pmSourceLocalSky(const psImage *image, const 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
+specified are used to measure the median \tbd{(clipped mean?)
+(statistic from 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
@@ -1242 +1250 @@
 
 \begin{verbatim}
-psSource *pmSourceMoments (psImage *image, psSource *source, float radius), 
+psSource *pmSourceMoments(psSource *source, const psImage *image, 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}
+\code{source.moments.sky} provided as the local background value.  The
+resulting moment values are applied to the \code{source.moments}
+entry, and the source is returned.  The moments are measured within
+the given radius of the \code{source.peak} coordinates.
+
+\begin{verbatim}
+psSource *pmSourceRoughClass(psSource *source, float saturate, float SNlim, const psRegion *valid);
+\end{verbatim}
+
+Based on the specified data values, make a guess at the source
+classification.  
+
+\subsection{Object Fitting}
+
+We need a way to fit a particular functional model to an object.
+PSLib includes the \code{psMinimizeLMChi2} and \code{psMinimizePowell}
+functions, which form the core of this processes.  However, additional
+support functions and wrapping functions are necessary for the
+specific case of source fitting.  The operations can be broken down
+into discrete steps:
+
+\begin{enumerate}
+\item Identify the pixels of interest
+
+\item Make a guess at the model parameters.  For some models, the
+parameters may be guessed based on only the moments.  For others,
+additional measurements must be made.
+
+\item Construct the input vectors from the pixels of interest.
+
+\item Apply fitting function \code{psMinimizeLMChi2()}
+
+\item Construct model image.
+
+\item Subtract model from image.
+\end{enumerate}
+
+\begin{verbatim}
+bool pmSourceSetPixelsCircle(psSource *source, const psImage *image, float radius);
+\end{verbatim}
+
+Define pixels associated with a source based on a circular aperture.
+This operation creates the \code{source.pixels} and \code{source.mask}
+entries for the source based on a circular aperture centered on the
+source centroid (or peak?).  The \code{source.pixels} is a subimage of
+the input image.  The function returns \code{TRUE} on success or
+\code{FALSE} on failure.
+
+\begin{verbatim}
+bool pmSourceModelGuess(psSource *source, const psImage *image);
+\end{verbatim}
+
+Convert available data to an initial guess for the given model.  The
+method of defining the model parameter guesses are specified for each
+model below.  The guess values are placed in the model parameters.  The
+function returns \code{TRUE} on success or \code{FALSE} on failure.
+
+\begin{verbatim}
+psArray *pmSourceContour(const psSource *source, const psImage *image, float level, int mode);
+\end{verbatim}
+
+Find points in a contour for the given source at the given level.  If
+mode is \code{PS_CONTOUR_CRUDE}, the contour is found by starting at
+the source peak, running along each pixel row until the level is
+crossed, then interpolating to the level coordinate for that row.
+This is done for each row, with the starting point determined by the
+midpoint of the previous row, until the starting point has a value
+below the contour level.  The resulting contour consists of two
+vectors giving the x and y coordinates of the contour levels.  This
+function may be used as part of the model guess inputs.
+
+\tbd{Other modes may be specified in the future for more refined contours}
+
+\begin{verbatim}
+bool pmSourceFitModel(psSource *source, psImage *image);
+\end{verbatim}
+
+Fit the requested model to the specified source.  The starting guess
+for the model is given by the input \code{source.model} parameter
+values.  The pixels of interest are specified by the
+\code{source.pixels} and \code{source.mask} entries.  This function
+calls \code{psMinimizeLMChi2()} on the image data.  The function
+returns \code{TRUE} on success or \code{FALSE} on failure.
+
+\begin{verbatim}
+bool pmSourceAddModel(psImage *image, psSource *source, bool center);
+\end{verbatim}
+
+Add the given source model flux to the provided image.  The boolean
+option center selects if the source is recentered to the image center
+or if it is placed at its centroid location.  The pixel range in the
+target image is at most the pixel range specified by the
+\code{source.pixels} image.  The success status is returned.  
+
+\begin{verbatim}
+bool pmSourceSubModel(psSource *source);
+\end{verbatim}
+
+Subtract the model from its image pixels given by
+\code{source.pixels}.  The success status is returned.  
 
 \subsection{Basic Object Models}
-
-float PGaussian (psVector *deriv, psVector *params, psVector *x);
-
-  X = x[0] - param[2];
-  Y = x[1] - param[3];
+\label{ObjectModels}
+
+We specify a variety of basic object models which are required.
+Details of the model functional forms, parameters, and the derivatives
+are specified in the ADD.
+
+\subsubsection{Real 2D Gaussian}
+
+\begin{verbatim}
+float psMinLM_Gauss2D(psVector *deriv, psVector *params, psVector *x);
+\end{verbatim}
+
+This function is a two-dimensional Gaussian with an elliptical
+cross-section and a constant local background.  
+
+The intial guess for the Gaussian parameters may be taken from the
+moments, peak value, and local sky.
+
+\subsubsection{Pseudo-Gaussian}
+
+\begin{verbatim}
+float psMinLM_PseudoGauss2D(psVector *deriv, psVector *params, psVector *x);
+\end{verbatim}
+
+This function is a polynomial approximation of a 2D Gaussian otherwise
+very similar to the real Gaussian.  It is used in place of a real
+Gaussian for speed.
+
+The intial guess for the Gaussian parameters may be taken from the
+moments, peak value, and local sky.
+
+\subsubsection{Waussian}
+
+\begin{verbatim}
+float psMinLM_Wauss2D(psVector *deriv, psVector *params, psVector *x);
+\end{verbatim}
+
+The Waussian is a modified polynomial approximation of a 2D Gaussian,
+with non-linear polynomial terms having variable coefficients, rather
+than the Taylor series values of 1/2 and 1/6.  
+
+\subsubsection{Twisted Gaussian}
+
+\begin{verbatim}
+float psMinLM_TwistGauss2D(psVector *deriv, psVector *params, psVector *x);
+\end{verbatim}
+
+This function describes an object with power-law wings and a flattened
+core, where the core has a different contour from the wings.  
+
+The intial guess for the Gaussian parameters may be taken from the
+moments, peak value, and local sky.
+
+\tbd{future galaxy models to be implemented}
+
+\subsubsection{Sersic Galaxy Model}
+
+\begin{verbatim}
+float psMinLM_Sersic(psVector *deriv, psVector *params, psVector *x);
+\end{verbatim}
+
+\subsubsection{Sersic with Core Galaxy Model}
+
+\begin{verbatim}
+float psMinLM_SersicCore(psVector *deriv, psVector *params, psVector *x);
+\end{verbatim}
+
+\subsubsection{Pseudo Sersic Galaxy Model}
+
+\begin{verbatim}
+float psMinLM_PseudoSersic(psVector *deriv, psVector *params, psVector *x);
+\end{verbatim}
+
+\appendix
+
+\section{Pseudo-C PSPhot}
+\label{psphot}
+
+\begin{verbatim}
+# include <pslib.h>
+# include <psmodule.h>
+
+main () {
+
+  psMetadata *header;
+  psImage *image;
+
+  fd = psFitsOpen (argv[1]);
+  md = psFitsReadHeader (fd);
+  image = psFitsReadImage (fd, md);
+
+  stats = psImageStats (NULL, image);
+
+  RDNOISE  = psMetadataLookup (md, "RDNOISE");
+  GAIN     = psMetadataLookup (md, "GAIN");
+  INNER    = psMetadataLookup (config, "INNER_RADIUS");
+  OUTER    = psMetadataLookup (config, "OUTER_RADIUS");
+  SATURATE = psMetadataLookup (config, "SATURATE");
+  NSIGMA   = psMetadataLookup (config, "PSF_PEAK_THRESHOLD");
+  RADIUS   = psMetadataLookup (config, "PSF_MOMENTS_RADIUS");
+  XBORDER  = psMetadataLookup (config, "XBORDER");
+  YBORDER  = psMetadataLookup (config, "YBORDER");
+
+  keep = psRegionAlloc (XBORDER, image->nCol - YBORDER, 
+                        YBORDER, image->nRow - YBORDER);
   
-  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];
+  Sky = stats->median;
+  Sig = sqrt(Sky/GAIN + SQ(RDNOISE));
+
+  kernel = psKernelParts ();
+  smooth = psImageConvolve (NULL, image, kernel, PS_PARTS);
+
+  peaks = pmFindImagePeaks (smooth, NSIGMA*Sig + Sky);
+
+  peaks = pmCullImagePeaks (peaks, SATURATE, keep);
   
-  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];
+  sources = pmSourceLocalSky (image, peaks, INNER, OUTER);
+
+  sources = pmSourceMoments (image, sources, RADIUS);
   
-  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;
-
-
+  sources = pmSourceRoughClassify (sources, SATURATE, MIN_SN_LIM, keep);
+
+  stars = pmSourceSelectBrightStars (sources);
+
+  stars = pmSourceFitModel (
+}
+
+
+
+psArray *pmFindImagePeaks (psImage *image, float threshold) {
+
+  psVector *row;
+
+  row = psVectorAlloc (image[0].Ncol);
+
+  /* find peaks in each row */
+  for (i = 0; i < image[0].Nrow; i++) {
+    rowpeaks = pmFindVectorPeaks (row, threshold);
+    peaks.x = rowpeaks;
+    peaks.y = i;
+    peaks.z = image (i, x);
+  }
+
+  /* drop non-local peaks (peaks with neighbors) */
+  for (n = 0; n < peaks.n; n++) {
+    if (!local_peak) {
+      drop_peak;
+    }
+  }
+  return (peaks);
+}
+\end{verbatim}
 
 \section{Revision Change Log}