Index: /trunk/doc/modules/ChangeLogSDRS.tex =================================================================== --- /trunk/doc/modules/ChangeLogSDRS.tex (revision 3069) +++ /trunk/doc/modules/ChangeLogSDRS.tex (revision 3070) @@ -1,3 +1,3 @@ -%%% $Id: ChangeLogSDRS.tex,v 1.10 2004-12-24 02:30:46 price Exp $ +%%% $Id: ChangeLogSDRS.tex,v 1.11 2005-01-22 01:59:10 eugene Exp $ \subsection{Changes from version 00 to version 01} @@ -35,3 +35,5 @@ \item Added a \code{mask} to \code{pmSubtractSky}, and specified that binned pixels which are clipped may be interpolated over, or simply ignored. \item Added further explanation for \code{pmReadoutCombine}. +\item Added Object Detection section +\item Added PSPhot pseudo-C example \end{itemize} Index: /trunk/doc/modules/ModulesSDRS.tex =================================================================== --- /trunk/doc/modules/ModulesSDRS.tex (revision 3069) +++ /trunk/doc/modules/ModulesSDRS.tex (revision 3070) @@ -1,3 +1,3 @@ -%%% $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,5 +11,5 @@ \project{Pan-STARRS Image Processing Pipeline} \organization{Institute for Astronomy} -\version{02} +\version{03} \docnumber{PSDC-430-012} @@ -27,6 +27,8 @@ % 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,5 +1042,5 @@ \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,5 +1049,5 @@ 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,65 +1055,66 @@ 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,29 +1123,16 @@ \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,27 +1142,29 @@ \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,5 +1177,5 @@ PS_MODEL_WAUSS; PS_MODEL_SERSIC; - PS_MODEL_SERSIC_BULGE; + PS_MODEL_SERSIC_CORE; } psModelType; \end{verbatim} @@ -1191,23 +1183,30 @@ \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,10 +1214,18 @@ 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,14 +1233,15 @@ 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,116 +1250,262 @@ \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 +# include + +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} Index: /trunk/doc/pslib/ChangeLogADD.tex =================================================================== --- /trunk/doc/pslib/ChangeLogADD.tex (revision 3069) +++ /trunk/doc/pslib/ChangeLogADD.tex (revision 3070) @@ -17,3 +17,4 @@ \item Added short section on histograms in the presence of errors. \item Added short note on inverse spherical transformations. +\item Added section on astronomical object models \end{itemize} Index: /trunk/doc/pslib/ChangeLogSDRS.tex =================================================================== --- /trunk/doc/pslib/ChangeLogSDRS.tex (revision 3069) +++ /trunk/doc/pslib/ChangeLogSDRS.tex (revision 3070) @@ -1,3 +1,3 @@ -%%% $Id: ChangeLogSDRS.tex,v 1.57 2005-01-20 08:25:45 price Exp $ +%%% $Id: ChangeLogSDRS.tex,v 1.58 2005-01-22 01:57:42 eugene Exp $ \subsection{Changes from version 00 to version 01} @@ -390,5 +390,5 @@ \end{itemize} -\subsection{Changes from Revision 10 (30 November 2004) to present} +\subsection{Changes from Revision 10 (30 November 2004) to Revision 11 (21 January 2005)} \begin{itemize} @@ -399,6 +399,5 @@ \item fixed error of psHash to psMetadata in psFitsReadHeaderSet \item added psFitsWriteImage -\item changed psFitsWriteImageSection to psFitsUpdateImage {\bf verify} -\item changed psFitsWriteHeader to psFitsUpdateHeader {\bf verify} +\item changed psFitsWriteImageSection to psFitsUpdateImage \item added header entry to psFitsWriteTable \item added psFitsUpdateTable @@ -434,3 +433,4 @@ PS_INTERPOLATE_BICUBIC_VARIANCE, PS_INTERPOLATE_SINC_VARIANCE}. \item Added \code{psImageTransform}. -\end{itemize} +\item Added section of Database Functions +\end{itemize} Index: /trunk/doc/pslib/psLibADD.tex =================================================================== --- /trunk/doc/pslib/psLibADD.tex (revision 3069) +++ /trunk/doc/pslib/psLibADD.tex (revision 3070) @@ -1,3 +1,3 @@ -%%% $Id: psLibADD.tex,v 1.56 2004-12-21 21:37:08 price Exp $ +%%% $Id: psLibADD.tex,v 1.57 2005-01-22 01:57:42 eugene Exp $ \documentclass[panstarrs]{panstarrs} @@ -1646,5 +1646,5 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsubsection{Missing and Todo} +\subsection{Missing and Todo} \tbd{define SINC, LAGRANGE interpolation} @@ -1659,4 +1659,263 @@ \tbd{define Brent's method \& minimization bracketing} + +\section{Pan-STARRS Modules} + +\subsection{Object Models} + +\subsubsection{Real 2D Gaussian} + +This function is a two-dimensional Gaussian with an elliptical +cross-section and a constant local background: +\[ +f(x,y) = Z_o e^{-z} + S_o +\] +where +\[ +z = \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y-y_o)^2}{2\sigma_y^2} + (x-x_o) (y - y_o) \sigma_{xy} +\] + +Below is the relationship between the \code{psModel} parameters and +the function parameters, sample C-code implementing the function +efficiently, and the value of the derivatives: + +\begin{verbatim} + param[0] = So; + param[1] = Zo; + param[2] = Xo; + param[3] = Yo; + param[4] = sqrt(2) / SigmaX; + param[5] = sqrt(2) / SigmaY; + param[6] = Sxy; + + 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 = exp(-z); + f = param[1]*r + param[0]; + /* f is the function value */ + + 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; +\end{verbatim} + +The intial guess for the Gaussian parameters may be taken from the +moments, peak value, and local sky. + +\subsubsection{Pseudo-Gaussian} + +This function is a polynomial approximation of a 2D Gaussian. The +function is very similar to the real Gaussian: +\[ +f(x,y) = Z_o (1 + z + z^2/2 + z^3/6)^{-1} + S_o +\] +where +\[ +z = \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y-y_o)^2}{2\sigma_y^2} + (x-x_o) (y - y_o) \sigma_{xy} +\] + +Below is the relationship between the \code{psModel} parameters and +the function parameters, sample C-code implementing the function +efficiently, and the value of the derivatives: + +\begin{verbatim} + param[0] = So; + param[1] = Zo; + param[2] = Xo; + param[3] = Yo; + param[4] = sqrt(2) / SigmaX; + param[5] = sqrt(2) / SigmaY; + param[6] = Sxy; + + 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; + t = 1 + z + 0.5*z*z; + r = 1.0 / (t*(1 + z/3)); /* ~ exp (-Z) */ + f = param[1]*r + param[0]; + /* f is the function value */ + + /* note difference from a pure gaussian: q = param[1]*r */ + q = param[1]*r*r*t; + 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; +\end{verbatim} + +The intial guess for the Gaussian parameters may be taken from the +moments, peak value, and local sky. + +\subsubsection{Waussian} + +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. The +function is very similar to the pseudo-Gaussian: +\[ +f(x,y) = Z_o (1 + z + B_2 (z^2/2 + B_3 z^3/6))^{-1} + S_o +\] +where +\[ +z = \frac{(x - x_o)^2}{2\sigma_x^2} + \frac{(y-y_o)^2}{2\sigma_y^2} + (x-x_o) (y - y_o) \sigma_{xy} +\] + +Below is the relationship between the \code{psModel} parameters and +the function parameters, sample C-code implementing the function +efficiently, and the value of the derivatives. Note the fudge factors +of 100 in the derivatives of $B_2$ and $B_3$: these are included to +slow the variation of these parameters, which are otherwise very +sensitive to small errors. + +\begin{verbatim} + param[0] = So; + param[1] = Zo; + param[2] = Xo; + param[3] = Yo; + param[4] = Sx; + param[5] = Sy; + param[6] = Sxy; + param[7] = B2; + param[8] = B3; + + 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; + t = 0.5*z*z*(1 + param[8]*z/3); + r = 1.0 / (1 + z + param[7]*t); /* ~ exp (-Z) */ + f = param[1]*r + param[0]; + + /* note difference from gaussian: q = param[1]*r */ + q = param[1]*r*r*(1 + param[7]*z*(1 + param[8]*z/2)); + 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*t; + deriv[8] = -100*param[1]*r*r*param[7]*(z*z*z)/6; + /* the values of 100 dampen the swing of param[7,8] */ +\end{verbatim} + +\subsubsection{Twisted Gaussian} + +This function describes an object with power-law wings and a flattened +core, where the core has a different contour from the wings. + +\[ +f(x,y) = Z_{\rm pk} (1 + z_1 + z_2^M)^{-1} + Sky +\] +where +\[ +z_1 = \frac{x^2}{2\sigma_{x,in}^2} + \frac{y^2}{2\sigma_{y,in}^2} + x y \sigma_{xy,in} +z_2 = \frac{x^2}{2\sigma_{x,out}^2} + \frac{y^2}{2\sigma_{y,out}^2} + x y \sigma_{xy,out} +\] + +\begin{verbatim} + param[0] = So; + param[1] = Zo; + param[2] = Xo; + param[3] = Yo; + param[4] = SxInner; + param[5] = SyInner; + param[6] = SxyInner; + param[7] = SxOuter; + param[8] = SyOuter; + param[9] = SxyOuter; + param[10] = N; + + 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,param[10])); + f = param[5]*r + param[6]; + + q1 = param[5]*SQ(r); + q2 = param[5]*SQ(r)*param[10]*pow(z2,(param[10]-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; + deriv[10] = -q1*ln(z2); +\end{verbatim} + +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} + +\begin{verbatim} +float Sersic() + param[0] = So; + param[1] = Zo; + param[2] = Xo; + param[3] = Yo; + param[4] = Sx; + param[5] = Sy; + param[6] = Sxy; + param[7] = Nexp; + +float SersicBulge() + param[0] So; + param[1] Zo; + param[2] Xo; + param[3] Yo; + param[4] SxInner; + param[5] SyInner; + param[6] SxyInner; + param[7] Zd; + param[8] SxOuter; + param[9] SyOuter; + param[10] = SxyOuter; + param[11] = Nexp; +\end{verbatim} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Index: /trunk/doc/pslib/psLibSDRS.tex =================================================================== --- /trunk/doc/pslib/psLibSDRS.tex (revision 3069) +++ /trunk/doc/pslib/psLibSDRS.tex (revision 3070) @@ -1,3 +1,3 @@ -%%% $Id: psLibSDRS.tex,v 1.168 2005-01-20 08:25:55 price Exp $ +%%% $Id: psLibSDRS.tex,v 1.169 2005-01-22 01:57:42 eugene Exp $ \documentclass[panstarrs,spec]{panstarrs} @@ -11,5 +11,5 @@ \project{Pan-STARRS Image Processing Pipeline} \organization{Institute for Astronomy} -\version{10} +\version{11} \docnumber{PSDC-430-007} @@ -41,4 +41,5 @@ 09 & 2004 Nov 15 & final for cycle 4 \\ \hline 10 & 2004 Nov 30 & update for cycle 4 \\ +11 & 2005 Jan 21 & draft for cycle 5 \\ \RevisionsEnd @@ -97,5 +98,5 @@ search path, PSLib may be used within a program by including the line \code{#include } into the C code and linking with -\code{-lpslib}. +\code{-lpslib.} This document describes the data structures and details the functions @@ -109,5 +110,5 @@ principle relevance and VerbPhase describes the operation applied to that data type. For example, the function which copies an image (of -type \code{psImage}) is called \code{psImageCopy()}. +type \code{psImage}) is called \code{psImageCopy().} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -4193,5 +4194,5 @@ typedef xmlDocPtr psXMLDoc; void psXMLDocFree(psXMLDoc *doc); -\end{verbatim}. +\end{verbatim} The next pair of functions convert a \code{psMetadata} data structure @@ -4200,5 +4201,5 @@ psXMLDoc *psMetadataToXMLDoc(const psMetadata *metadata); psMetadata *psXMLDocToMetadata(const psXMLDoc *doc); -\end{verbatim}. +\end{verbatim} The next pair of functions loads the data in a named file into a @@ -4208,5 +4209,5 @@ psXMLDoc *psXMLParseFile(const char *filename); int psXMLDocToFile(const psXMLDoc *doc, const char *filename); -\end{verbatim}. +\end{verbatim} The next pair of functions accepts a block of memory and parses it @@ -4215,5 +4216,5 @@ psXMLDoc *psXMLParseMemory(const char *buffer, const int size); int psXMLDocToMemory(const psXMLDoc *doc, char *buffer); -\end{verbatim}. +\end{verbatim} The next pair of functions read from and write to a file descriptor. @@ -4224,5 +4225,174 @@ psXMLDoc *psXMLParseFD(int fd); int psXMLDocToFD(const psXMLDoc *doc, int fd); -\end{verbatim}. +\end{verbatim} + +\subsection{Database Functions} + +Many of the applications that PSLib will be used for will require +access to a simple relational database. PSLib includes generic +database-independent interface mechanisms as part of its API set. The +most important aspect of PSLib's database support is to abstract as +much database specific complexity as is feasible. As almost all RDBMS +provide at least a simple transactional model, commit and rollback +support should be provided. + +Currently, only support for MySQL 4.1.x is required but other backends +may be added as options in the future. As a particular example which +has implications for the database interaction model, support for +SQLite may be required in the future. Currently, the choice of +backend database interface may be made as a compile option. Details +of the specified APIs in the discussion below refer to the relevant +MySQL functions. + +Database errors must be trapped and placed onto the psError stack. +The complete error message should be retrieved with the database's +error function. + +\subsubsection{Managing the Database Connection} + +We specify a database handle which carries the information about the +database connection: + +\begin{verbatim} + typedef struct { + MYSQL mysql; + } psDB; +\end{verbatim} + +The following collection of functions provides basic database functionality: + +\begin{verbatim} + // wraps mysql_init() & mysql_real_connect() + psDB *psDBInit(const char *host, const char *user, const char *passwd, const char *dbname); + + // wraps mysql_close() + void psDBCleanup(psDB *dbh); + + // wraps mysql_create_db() + bool psDBCreate(psDB *dbh, const char *dbname); + + // wraps mysql_select_db() + bool psDBChange(psDB *dbh, const char *dbname); + + // wraps mysql_drop_db() + bool psDBDrop(psDB *dbh, const char *dbname); +\end{verbatim} + +For MySQL support, \code{psDBInit()} wraps \code{mysql_init()} and +\code{mysql_real_connect()} in order to initialize a psDB structure and +establish a database connection. A null pointer should be returned on +failure. + +When implementing support for SQLite, or other DB which is purely +file-based, the \code{host}, \code{user}, and \code{passwd} arguments +would be ignored while \code{dbname} would specify the path to the +SQLite db file. + +\subsubsection{Interacting with Database Tables} + +The functions in this section perform high level interactions with the +database tables. All of them should behave ``atomically'' with +respect to the state of the database. Specifically, all interactions +with the database should be done as a part of a transaction that is +rolled-back on failure and committed only after all queries used by +the API have been run. In general, this API set attempts to treat a +database table as a 2D matrix where columns can be represented by a +\code{psVector} and rows as a \code{psMetadata} type. A +\code{psMetadata} collection is also used to define the columns of a +table and as part of the query restrictions. + +\begin{verbatim} + bool psDBCreateTable(psDB *dbh, const char *tableName, psMetadata *md); +\end{verbatim} + +This function generates and executes the SQL needed to create a table +named \code{tableName}, with the column names and datatypes as +described in \code{md}. Each data item in the \code{psMetadata} +collection represents a single table field. The name of the field is +given by the name of the \code{psMetadataItem} and the data type is +give by the \code{psMetadataItem.type} and \code{psMetadataItem.ptype} +entries. A lookup table should be used to convert from PSLib types +into MySQL compatible SQL data types. For example, a +\code{PS_META_STR} would map to an SQL99 varchar. If value of +\code{type} is \code{PS_META_STR} then the \code{psMetadataItem.data} +element is set to a string with the length for the field written as a +text string. The value of the \code{psMetadataItem.data} element is +unused for the \code{PS_META_PRIMITIVE} types. Other metadata types +beyond \code{PS_META_STR} and \code{PS_META_PRIMITIVE} are not allowed +in a table definition metadata collection. + +Database indexes can be specified setting the \code{comment} field to +``\code{Primary Key}'' or ``\code{Key}''. Comment are otherwise +ignored. + +\begin{verbatim} + bool psDBDropTable(psDB *dbh, const char *tableName); +\end{verbatim} + +This function deletes the specified table. + +\begin{verbatim} + psArray *psDBSelectColumn(psDB *dbh, const char *tableName, const char *col, const psU64 limit); + psVector *psDBSelectColumnNum(psDB *dbh, const char *tableName, const char *col, const psU64 limit); +\end{verbatim} + +These functions generates and executes the SQL needed to select an +entire column from a table or up to \code{limit} rows from it. If +\code{limit} is 0, the entire range is returned. The database +response is processed and a \code{psArray} of elements of the native +type is returned. The Num version of the function returns the data in +a \code{psVector}. It returns an error (NULL) if the requested field +is not a numerical type. + +\begin{verbatim} + psArray *psDBSelectRows(psDB *dbh, const char *tableName, psMetadata *where); +\end{verbatim} + +This function returns rows from the specified table which match +the restrictions given by \code{where}. The restrictions are +specified as field / value pairs. The \code{psMetadata} collection +where must consist of valid database fields, though the database query +checking functions may be used to validate the fields as part of the +query. If \code{where} is \code{NULL}, then there are no restrictions +on the rows selected. The selected rows are returned as a +\code{psArray} of \code{psMetadata} values, one per row. + +\begin{verbatim} + bool psDBInsertRow(psDB *dbh, const char *tableName, psMetadata *row); +\end{verbatim} + +Insert the data from \code{row} into \code{tableName}. It should be noted in +the API reference that if fields are specified in \code{row} that do not exist +in \code{tablename}, the insert will fail. + +\begin{verbatim} + psArray *psDBDumpRows(psDB *dbh, const char *tableName); +\end{verbatim} + +Fetch all rows as an psArray of psMetadata. + +\begin{verbatim} + psArray *psDBDumpCols(psDB *dbh, const char *tableName); +\end{verbatim} + +Fetch all columns as an psArray of psVector + +\begin{verbatim} + psU64 psDBUpdateRow(psDB *dbh, const char *tableName, psMetadata *where, psMetadata *values); +\end{verbatim} + +Update the columns contained in \code{values} in the row(s) that have a field +with the value indicated by \code{where} (note that this is only allows very +limited use of SQL99's ``where'' semantics). The number of rows modified is +returned. If there are multiple psMetadataItems in \code{where} then each item +should be considered as an additional constraint. e.g. ``where foo = x and +where bar = y'' + +\begin{verbatim} + psU64 psDBDeleteRow(psDB *dbh, const char *tableName, psMetadata *where); +\end{verbatim} + +Delete the rows that are matched by \code{where} using the same semantics for +\code{where} as in psDBUpdateRow(). \subsection{FITS I/O Functions} @@ -4598,10 +4768,7 @@ \begin{verbatim} -psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out, const psPlaneTransform *in, float xMin, - float xMax, float yMin, float yMax, int nSamples); -psPlaneTransform *psPlaneTransformCombine(psPlaneTransform *out, const psPlaneTransform *trans1, - const psPlaneTransform *trans2); -bool psPlaneTranformFit(psPlaneTransform *trans, const psArray *source, const psArray *dest, int nRejIter, - float sigmaClip); +psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out, const psPlaneTransform *in, psRegion *region, int nSamples); +psPlaneTransform *psPlaneTransformCombine(psPlaneTransform *out, const psPlaneTransform *trans1, const psPlaneTransform *trans2); +bool psPlaneTranformFit(psPlaneTransform *trans, const psArray *source, const psArray *dest, int nRejIter, float sigmaClip); \end{verbatim} @@ -4612,24 +4779,22 @@ the forward transformation. In the event that the input transformation is linear, an exact solution may be calculated; -otherwise \code{nSamples} samples in each axis, ranging from -\code{xMin} to \code{xMax} and \code{yMin} to \code{yMax} shall be -used as a grid to fit the best inverse transformation. The function -shall return \code{NULL} if it was unable to generate the inverse -transformation; otherwise it shall return the inverse transformation. -In the event that \code{out} is \code{NULL}, a new -\code{psPlaneTransform} shall be allocated and returned. - -\code{psPlaneTransformSubsume} takes two transformations +otherwise \code{nSamples} samples in each axis, covering the region +specified by \code{region} shall be used as a grid to fit the best +inverse transformation. The function shall return \code{NULL} if it +was unable to generate the inverse transformation; otherwise it shall +return the inverse transformation. In the event that \code{out} is +\code{NULL}, a new \code{psPlaneTransform} shall be allocated and +returned. + +\code{psPlaneTransformCombine} takes two transformations (\code{trans1} and \code{trans2}) and returns a single transformation that has the effect of performing \code{trans1} followed by \code{trans2}. The function shall return \code{NULL} if it was unable to generate the transformation; otherwise it shall return the -transformation. \tbd{Not sure on the algorithm yet --- it may be the -same as for \code{psPlaneTransformInvert}, in which case we will need -the ranges and number of samples as well.} - -\code{psPlaneTransform} takes two arrays containing matched -coordinates (i.e., coordinates the \code{source} correspond to the -appropriate coordinates in the \code{dest}) and returns the +transformation. + +\code{psPlaneTransformFit} takes two arrays containing matched +coordinates (i.e., coordinates in the \code{source} array correspond +to the coordinates in the \code{dest} array) and returns the best-fitting transformation. The \code{source} and \code{dest} will contain \code{psCoord}s. In the event that the number of coordinates