%%% $Id: ModulesSDRS.tex,v 1.81 2007-02-03 00:12:36 eugene Exp $
\documentclass[panstarrs,spec]{panstarrs}

% basic document variables
\title{Pan-STARRS PS-1 Image Processing Pipeline Modules}
\subtitle{Supplementary Design Requirements}
\shorttitle{Modules SDRS}
\author{Paul Price, Eugene Magnier}
\audience{Pan-STARRS PMO}
\group{Pan-STARRS Algorithm Group}
\project{Pan-STARRS Image Processing Pipeline}
\organization{Institute for Astronomy}
\version{11}
\docnumber{PSDC-430-012}

\begin{document}
\maketitle
\sloppy

% -- Revision History --
% provide explicit values for the old versions
% use '\theversion' for the current version (set above)
% use \hline between each table row
\RevisionsStart
% version  Date            Description
DR & 2004 Jun 7 & Draft \\ \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
04 & 2005 Feb 14 & final for cycle 5 \\ \hline
05 & 2005 Mar 21 & draft for cycle 6 \\ \hline
06 & 2005 Apr 27 & final for cycle 6 \\ \hline
07 & 2005 Jul 15 & final for cycle 7 \\ \hline
08 & 2005 Sep 13 & final for cycle 8 \\ \hline
09 & 2005 Oct 18 & final for cycle 9 \\ \hline
10 & 2005 Dec 06 & draft for cycle 10 \\ \hline
11 & 2006 Jan 22 & revisions for CDR \\ \hline
\RevisionsEnd

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\DocumentsInternal
PSCD-230-001  &   PS-1 Design Reference Mission \\ \hline
PSDC-430-004  &   Pan-STARRS PS-1 IPP C Code Conventions \\ \hline
PSDC-430-005  &   Pan-STARRS PS-1 IPP Software Requirements Specification \\ \hline
PSDC-430-006  &   Pan-STARRS PS-1 IPP Algorithm Design Document \\ \hline
PSDC-430-011  &   Pan-STARRS PS-1 IPP System/Subsystem Design Description \\ \hline
PSDC-430-???  &   Pan-STARRS PS-1 IPP Configurations \\ \hline
\DocumentsExternal
Posix Standard & Open Group Based Specifications Issue 6, IEEE Std 1003.1, 2003 \\
\DocumentsEnd

\tableofcontents
\pagebreak 

\listoffigures
\pagebreak 
\pagenumbering{arabic}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

This document describes the Pan-STARRS Image Processing Pipeline (IPP)
data analysis Modules.  The Modules use the functionality of the
Pan-STARRS Library (PSLib) to perform more complex tasks, especially
tasks which require assumptions of astronomical analysis or the data
organization.  Within the IPP, the Modules are tied together into
programs which perform complete data analysis tasks (an ``analysis
stage'').  The modules may be tied together within a C framework or
using a high-level scripting language.

In order to preserve name space, globally-visible structures and
functions shall be prefixed with \code{pm}, for ``Pan-STARRS
Modules''.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Configuration}

Configuration information in the IPP is provided on four levels: the
site, camera, format and recipe configurations.  The use of these, and
their contents are described in an additional document, PSDC-430-???,
``Pan-STARRS PS1 IPP Configurations''.  Here we define functions that
read these configuration files.


\subsection{Configuration information}

This structure stores the configuration information: the site, camera
and recipe configuration, the command-line arguments, the pmFPAfiles
used, and the database handle.
\begin{datatype}
typedef struct {
    psMetadata *site;                   ///< Site configuration
    psMetadata *camera;                 ///< Camera specification
    const char *cameraName;             ///< Camera name
    psMetadata *recipes;                ///< Recipes for processing
    psMetadata *arguments;              ///< Processed command-line arguments
    psMetadata *files;                  ///< pmFPAfiles used for analysis
    psDB *database;                     ///< Database handle
    int *argc;                          ///< Number of command-line arguments
    char **argv;                        ///< Command-line arguments (raw version)
    // Private members
    p_pmRecipeSource recipesRead;       ///< Which recipe sources have been read
    psMetadata *recipesSource;          ///< Where each recipe came from
} pmConfig;
\end{datatype}

The site, camera and recipe configurations are represented here.  The
camera format configuration is stored elsewhere (within the
\code{pmFPA}).

Note that the \code{arguments} are different from the \code{argc} and
\code{argv}, since it is the parsed version.  This metadata may also
be used to carry program options, if desired.

The allocator for this structure is:
\begin{prototype}
pmConfig *pmConfigAlloc(int *argc,      /// Number of command-line arguments
                        char **argv     /// Command-line arguments
                       );
\end{prototype}


\subsection{Reading configuration files}

Read configuration information from the command line:
\begin{prototype}
pmConfig *pmConfigRead(int *argc,       ///< Number of command-line arguments
                       char **argv      ///< Array of command-line arguments
                      );
\end{prototype}
pmConfigRead loads the site configuration (the file name is specified
by "-site SITE\_FILE" on the command-line, the PS\_SITE environment variable, or
it is \code{$HOME/.ipprc}).  The configuration search %$ path is
set. The camera configuration is loaded if it is specified on the
command line ("-camera CAMERA\_FILE"). Recipes specified on the command
line ("-recipe RECIPE\_NAME RECIPE\_SOURCE") are also loaded.  These
command-line arguments are removed from from the command-line, to
simplify parsing.  The psLib log, trace and time setups are also
performed if specified in the site configuration.

Read a configuration file:
\begin{prototype}
bool pmConfigFileRead(psMetadata **config, ///< Config to output
                      const char *name, ///< Name of file
                      const char *description ///< Description of file
                     );
\end{prototype}
Read a metadata configuration file into the supplied metadata.
Produce an error and return false if there's a problem.

Set static configuration information:
\begin{prototype}
void pmConfigSet(const char *path ///< Search paths for configuration files; colon-delimited directories
                );
\end{prototype}
The search path for the configuration files is a local static
variable, set by this function.

Free static memory used in the configuration system:
\begin{prototype}
void pmConfigDone(void);
\end{prototype}


\subsection{Camera configuration}

Validate a header against the camera format:
\begin{prototype}
bool pmConfigValidateCameraFormat(const psMetadata *cameraFormat, ///< Camera format containing the RULE
                                  const psMetadata *header // FITS header for the PHU
                                 );
\end{prototype}
Given a FITS header (the PHU header), check it against the RULE
metadata contained within the camera format; return true if it
matches.

Determine the camera format (and camera if unknown) from examining the
header:
\begin{prototype}
psMetadata *pmConfigCameraFormatFromHeader(pmConfig *config, ///< The configuration
        const psMetadata *header ///< The FITS header
                                          );
\end{prototype}
Given a FITS header, check it against all known cameras (unless we
already know which camera, from pmConfigRead) and all known formats
for those cameras in order to identify which is appropriate.  The
first matching format is accepted; further matches produce warnings.
The accepted camera is saved in the configuration.  The accepted
format is returned.

Return the camera configuration specified by name:
\begin{prototype}
psMetadata *pmConfigCameraByName(pmConfig *config, ///< The configuration
                                 const char *cameraName ///< The camera name header
                                );
\end{prototype}
Given a camera name, returns the camera configuration metadata.

Make the supplied header conform to the nominated camera format:
\begin{prototype}
bool pmConfigConformHeader(psMetadata *header, ///< Header to conform
                           const psMetadata *format ///< Camera format
                          );
\end{prototype}
Given a FITS header, make it conform to the RULE in the specified
camera format.  This is useful for switching between formats, or
generating fake data that must be recognised by
pmConfigCameraFormatFromHeader.


\subsection{Recipes}

Recipes are read from a variety of sources: the command line, the site
configuration, and the camera configuration.  In addition, recipes
specified on the command line may be a symbolic link to recipes
defined in the other locations.  Since the various sources of recipes
may be read at different times (e.g., the camera configuration may be
read from the command line, or it may only be read once a FITS file
has been opened for inspection), we specify the following list of
recipe sources:
\begin{datatype}
typedef enum {
    PM_RECIPE_SOURCE_NONE        = 0x00, ///< None yet
    PM_RECIPE_SOURCE_SITE        = 0x01, ///< Site configuration
    PM_RECIPE_SOURCE_CAMERA      = 0x02, ///< Camera configuration
    PM_RECIPE_SOURCE_CL          = 0x04, ///< Command-line
    PM_RECIPE_SOURCE_SYMBOLIC    = 0x14, ///< Symbolic link, specified on command-line
    PM_RECIPE_SOURCE_ALL         = 0xff  ///< All sources
} pmRecipeSource;
\end{datatype}

See PSDC-430-???, ``IPP Configurations'' for information on symbolic
links, and how these are identified and resolved, and also for the
precedence of sources.

Read recipes:
\begin{prototype}
bool pmConfigReadRecipes(pmConfig *config ///< Configuration
                        );
\end{prototype}
Attempt to read recipes from the sources that are available but have
not already been read.Having read a recipe, attempt to resolve
symbolic links that were specified on the command line.


\subsection{Database}

Setup the database:
\begin{prototype}
psDB *pmConfigDB(pmConfig *config       ///< Configuration
                );
\end{prototype}
Initialise the database connection using the DBSERVER, DBNAME, DBUSER,
DBPASSWORD values provided in the site configuration.  Stores the
database handle in the configuration, and also returns it.


\subsection{File sets}

Read the command-line for files (or a text file containing a list of
files):
\begin{prototype}
psArray *pmConfigFileSets(pmConfig *config, ///< Configuration, containing command-line arguments
                          const char *file, ///< CL argument specifying a filename
                          const char *list ///< CL argument specifying a text file with a list of filenames
                         );
\end{prototype}
Given the 'file' and 'list' arguments (e.g., "-file" and "-list"),
find the arguments associated with these words and interpret them as
lists of files.  Return an array of the resulting filenames.

Stuff associated files from the command-line into a metadata:
\begin{prototype}
bool pmConfigFileSetsMD(psMetadata *metadata, ///< Metadata into which to stuff the array
                        pmConfig *config, ///< Configuration (which command-line arguments)
                        const char *name, ///< Name for array in the metadata
                        const char *file, ///< CL argument specifying a filename
                        const char *list ///< CL argument specifying a text file with a list of filenames
                       );
\end{prototype}
Calls pmConfigFileSets to parse the command line for filenames (or a
list which provides filenames), and stuffs the array of filenames into
the metadata under "name".

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{``Concepts''}
\label{sec:concepts}

Each image from an astronomical instrument has associated with it what
we will call {\it concepts} (for want of a better word; \tbd{we would
like to call this ``metadata'', but unfortunately that name is already
taken}).  These are values corresponding to general quantities and
qualities necessary to understand and interpret the data, such as
airmass, date, read noise and filter.  The values of each of the below
concepts shall be determined when the FPA is read into memory (via
\code{pmFPARead}), and stored at the appropriate level in the focal
plane hierarchy.

After ingest, the user may safely assume that all of the above
concepts exist at the appropriate level (meaning the user needn't be
hampered by excessive error checking), is of the specified type
(meaning the user doesn't need to worry about whether the value of
interest is stored in, e.g., floating point or double precision or
even a colon-delimited string) and in the specified format (meaning
the user doesn't need to know, e.g., whether the right ascension is in
radians or degrees) --- all the conversions are handled by the
``concepts'' functions at ingest.

Most of the structures and functions in this section are intended to
be ``private'', since there is no need envisioned for the user to call
them directly.

\subsection{Specifying a concept}

Specifying a ``concept'' requires a (meaningful) name (preferably with
the level in the name, e.g., \code{CELL.EXPOSURE}), a
comment/description, a type, a default or blank value, functions to
read and write, and a level that the concept applies to
(FPA/Chip/Cell).

\begin{datatype}
typedef struct {
    psMetadataItem *blank;              ///< Blank value of concept; also contains the name and comment
    pmConceptParseFunc parse;           ///< Function to call to read the concept
    pmConceptFormatFunc format;         ///< Function to call to write the concept
} pmConceptSpec;
\end{datatype}

The allocator is:
\begin{prototype}
pmConceptSpec *pmConceptSpecAlloc(psMetadataItem *blank, ///< Blank value; contains the name
                                  pmConceptParseFunc parse, ///< Function to call to parse the concept
                                  pmConceptFormatFunc format ///< Function to call to format the concept
                                 );
\end{prototype}

\code{blank} provides the name, type and default/blank value for the
concept.  \code{parse} and \code{format} provide the functions to
parse and format the concepts:

\begin{datatype}
typedef psMetadataItem* (*pmConceptParseFunc)(const psMetadataItem *concept, ///< Concept to parse
        const psMetadataItem *pattern, ///< Pattern for parsing
        const psMetadata *cameraFormat, ///< Camera format definition
        const pmFPA *fpa, ///< FPA for concept, or NULL
        const pmChip *chip, ///< Chip for concept, or NULL
        const pmCell *cell ///< Cell for concept, or NULL
                                             );
typedef psMetadataItem* (*pmConceptFormatFunc)(const psMetadataItem *concept, ///< Concept to format
        const psMetadata *cameraFormat, ///< Camera format definition
        const pmFPA *fpa, ///< FPA for concept, or NULL
        const pmChip *chip, ///< Chip for concept, or NULL
        const pmCell *cell ///< Cell for concept, or NULL
                                              );
\end{datatype}

A concept that has just been read requires parsing, so that the final
product stored in the \code{concepts} metadata is of the correct type
and format.  The parsing function requires a \code{pattern} (generally
the \code{blank}), which provides the name and comment for the parsed
concept.  The camera format (\code{cameraFormat}) allows looking up of
the \code{FORMATS} or other information, and the \code{fpa},
\code{chip} and \code{cell} allow lookup of other concepts.

Before writing, concepts are formatted so that the source's format
(likely differing from that of the psModules concepts) may be
respected.  Again, providing the camera format (\code{cameraFormat}),
\code{fpa}, \code{chip} and \code{cell} allow the formatting function
to access information necessary to format.


\subsection{Registering a concept}

The concept specifications that have been registered are stored in
three (static) \code{psMetadata}s, one for each level (FPA, chip,
cell).

Registering a concept is achieved by:
\begin{prototype}
bool pmConceptRegister(psMetadataItem *blank, ///< Blank value; contains the name and default comment
                       pmConceptParseFunc parse, ///< Function to call to parse the concept, or NULL
                       pmConceptFormatFunc format, ///< Function to call to format the concept, or NULL
                       pmFPALevel level ///< Level at which to store concept in the FPA hierarchy
                      );
\end{prototype}

\code{pmConceptRegister} generates a concept specification from the
provided \code{blank}, and \code{parse} and \code{format} functions,
and register it in the metadata specified by the \code{level}.

A \code{parse} function of \code{NULL} indicates that there is no
special interpretation of the concept required, and that it can be
used as read.  A \code{write} function of \code{NULL} indicates that
no special formatting of the concept is required, and that it can be
written as is.


\subsection{Default concepts}

Below is a list of concepts that are defined in the psModules, with
the expected type and a short description.

\begin{itemize}
\item \code{FPA.CAMERA} (STR): Camera used
\item \code{FPA.FOCUS} (F32): Telescope focus
\item \code{FPA.AIRMASS} (F32): Airmass at boresight
\item \code{FPA.FILTER} (STR): Filter used
\item \code{FPA.POSANGLE} (F32): Position angle of instrument
\item \code{FPA.RADECSYS} (STR): Celestial coordinate system
\item \code{FPA.RA} (F64): Right Ascension of boresight
\item \code{FPA.DEC} (F64): Declination of boresight
\item \code{FPA.OBSTYPE} (STR): Type of observation
\item \code{FPA.OBJECT} (STR): Object of observation
\item \code{FPA.ALT} (F64): Altitude of telescope
\item \code{FPA.AZ} (F64): Azimuth of telescope
\item \code{FPA.TIMESYS} (STR): Time system
\item \code{FPA.TIME} (psTime): Time of exposure
\item \code{CHIP.XPARITY} (S32): Orientation in x compared to the rest of the FPA, -1 or +1
\item \code{CHIP.YPARITY} (S32): Orientation in y compared to the rest of the FPA, -1 or +1
\item \code{CHIP.X0} (S32): Position of (0,0) on the FPA
\item \code{CHIP.Y0} (S32): Position of (0,0) on the FPA
\item \code{CHIP.TEMP} (F32): Temperature of chip
\item \code{CELL.GAIN} (F32): CCD gain (e/count)
\item \code{CELL.READNOISE} (F32): CCD read noise (e)
\item \code{CELL.SATURATION} (F32): Saturation level (counts)
\item \code{CELL.BAD} (F32): Bad level (counts)
\item \code{CELL.XPARITY} (S32): Orientation in x compared to the rest of the chip, -1 or +1
\item \code{CELL.YPARITY} (S32): Orientation in y compared to the rest of the chip, -1 or +1
\item \code{CELL.READDIR} (S32): Read direction, rows=1, cols=2
\item \code{CELL.EXPOSURE} (F32): Exposure time (sec)
\item \code{CELL.DARKTIME} (F32): Time since flush (sec)
\item \code{CELL.TRIMSEC} (psRegion): Trim section
\item \code{CELL.BIASSEC} (psList of psRegions): Bias sections
\item \code{CELL.XBIN} (S32): Binning in x
\item \code{CELL.YBIN} (S32): Binning in y
\item \code{CELL.TIMESYS} (STR): Time system
\item \code{CELL.TIME} (psTime): Time of exposure
\item \code{CELL.X0} (S32): Position of (0,0) on the chip
\item \code{CELL.Y0} (S32): Position of (0,0) on the chip
\end{itemize}

\tbd{CELL.EXPOSURE, CELL.DARKTIME and CELL.TIME should actually be
specified at the readout level.  However, at this present time, we're
not sure how these should be specified, and so we move them up to the
cell level and assume that all readouts are of the same exposure and
dark time.}

The concept specifications for the above may be registered by
\code{pmConceptsInit}:
\begin{prototype}
bool pmConceptsInit(void);
\end{prototype}

Since defined concept specifications are required before any concept
ingest can take place, all functions that work with the concepts must
call \code{pmConceptsInit} first.  This is usually performed internal
to the pmConcepts, so the user need not call this function.

The concept specification metadata containers and the concept
specifications that have been registered are freed by:
\begin{prototype}
void pmConceptsDone(void);
\end{prototype}
Calling \code{pmConceptsDone} is required in order to avoid a memory
leak, since the metadata containers are defined \code{static}.


\subsection{Concept sources}

Concepts may be read from, and written to a variety of sources: the
\code{CELLS} information within the camera format, the FITS header
(either the primary header, or extension), default values from the
camera format, or the database.

Since some sources become available at different times from others, we
need to provide some specificity to reading and writing concepts (or
we're forced to wait until everything's available, which we don't want
to do).  Concepts may be read from or written to multiple sources at
once by OR-ing them.

\begin{datatype}
typedef enum {
    PM_CONCEPT_SOURCE_NONE     = 0x00,  ///< No concepts
    PM_CONCEPT_SOURCE_CELLS    = 0x01,  ///< Concept comes from the camera information
    PM_CONCEPT_SOURCE_DEFAULTS = 0x02,  ///< Concept comes from defaults
    PM_CONCEPT_SOURCE_PHU      = 0x04,  ///< Concept comes from PHU
    PM_CONCEPT_SOURCE_HEADER   = 0x08,  ///< Concept comes from FITS header
    PM_CONCEPT_SOURCE_DATABASE = 0x10,  ///< Concept comes from database
    PM_CONCEPT_SOURCE_ALL      = 0xff   ///< All concepts
} pmConceptSource;
\end{datatype}


\subsection{Blanking}

By ``blanking'', we mean setting the concepts to a default or blank
value (e.g., \code{NaN} for floating point); this takes place before
reading, and can be used to set up a focal plane hierarchy without
reading from any particular source.  This is not normally done by the
user, but it is performed on the focal plane hierarchy at construction
(see \code{pmFPAConstruct}).

The following functions blank the concepts at the appropriate level:

\begin{prototype}
bool pmConceptsBlankFPA(pmFPA *fpa);
bool pmConceptsBlankChip(pmChip *chip);
bool pmConceptsBlankCell(pmCell *cell);
\end{prototype}


\subsection{Reading}

Reading a concept is the act of looking up the value, parsing it, and
setting the parsed concept in the \code{concepts} metadata in the
focal plane hierarchy at the appropriate level.  It is not normally
necessary for the user to call these functions, since they are called
when constructing or reading the FPA (see \code{pmFPAConstruct},
\code{pmFPARead}).

The following functions read the concepts at the appropriate level:

\begin{prototype}
bool pmConceptsReadFPA(pmFPA *fpa,      ///< FPA for which to read concepts
                       pmConceptSource source, ///< Source for concepts
                       bool propagateDown, ///< Propagate to lower levels?
                       psDB *db         ///< Database handle
                      );
bool pmConceptsReadChip(pmChip *chip,   ///< Chip for which to read concepts
                        pmConceptSource source, ///< Source for concepts
                        bool propagateUp, ///< Propagate to higher levels?
                        bool propagateDown, ///< Propagate to lower levels?
                        psDB *db        ///< Database handle
                       );
bool pmConceptsReadCell(pmCell *cell,   ///< Cell for which to read concepts
                        pmConceptSource source, ///< Source for concepts
                        bool propagateUp, ///< Propagate to higher levels?
                        psDB *db        ///< Database handle
                       );

\end{prototype}
In each case, we require the hierarchy component for which to read the
concepts, an indication of the source from which to read concepts
(which source has recently become available), the database handle (for
reading concepts from the database), and flags indicating whether we
should propagate the read to higher or lower levels in the hierarchy.

An additional function reads concepts for as many hierarchy components
as are supplied:
\begin{prototype}
bool pmConceptsRead(pmFPA *fpa,         ///< FPA for which to read concepts
                    pmChip *chip,       ///< Chip for which to read concepts, or NULL
                    pmCell *cell,       ///< Cell for which to read concepts, or NULL
                    pmConceptSource source, ///< The source of the concepts to read
                    psDB *db            ///< Database handle
                   );
\end{prototype}


\subsection{Writing}

Writing a concept is the act of formatting the concept for the native
format, and writing it to the appropriate place (here, ``writing''
often means preparing for writing, by updating the header metadata, or
checking the concept value against the default).  The user will not
normally use these functions directly, since they are called when
writing the FPA (see \code{pmFPAWrite}).

\begin{prototype}
bool pmConceptsWriteFPA(const pmFPA *fpa,     ///< FPA for which to write concepts
                        pmConceptSource source, ///< Source for concepts
                        bool propagateDown, ///< Propagate to lower levels?
                        psDB *db        ///< Database handle
                       );
bool pmConceptsWriteChip(const pmChip *chip,  ///< Chip for which to write concepts
                         pmConceptSource source, ///< Source for concepts
                         bool propagateUp,///< Propagate to higher levels?
                         bool propagateDown, ///< Propagate to lower levels?
                         psDB *db       ///< Database handle
                        );
bool pmConceptsWriteCell(const pmCell *cell,  ///< FPA for which to write concepts
                         pmConceptSource source, ///< Source for concepts
                         bool propagateUp, ///< Propagate to higher levels?
                         psDB *db       ///< Database handle
                        );
\end{prototype}
In each case, we require the hierarchy component for which to read the
concepts, an indication of the source to which to write concepts
(which source requires updating), the database handle (for writing
concepts to the database), and flags indicating whether we should
propagate the write to higher or lower levels in the hierarchy.


\subsection{Copying concepts}

The values of concepts may be copied from one source to another:
\begin{prototype}
bool pmFPACopyConcepts(pmFPA *target, const pmFPA *source);
\end{prototype}

\code{pmFPACopyConcepts} iterates through the focal plane hierarchy,
copying the values of the concepts from the \code{source} to the
\code{target}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\input{CameraImages.tex}

%\input{CameraGeometry.tex}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Image Detrending}

Image Detrending is the image analysis process wherein the
instrumental signatures are removed from the individual images.  This
section discusses the modules used for image detrending.  The basic
image detrending steps are:
\begin{itemize}
\item Subtract overscan, bias and dark;
\item Correct for non-linearity;
\item Shutter correction;
\item Flat-field;
\item Mask bad pixels;
\item Fringe subtraction;
\item Subtract the background;
\item Mask cosmic rays;
\item Mask optical defects;
\end{itemize}

We also include here some of the (more complicated) functions for
generating the image detrend calibrations (fringes, flat-field
normalisation), and determining the suitable calibrations from the
detrend database.

\subsection{Bias subtraction}
\label{sec:bias}

The bias subtraction module provides a facility to correct detector
images for the electronic pedestal introduced by the readout
electronics.  The pedestal is measured from one or more of the
overscan, bias frame and dark frames.

The following function subtracts the overscan, bias and/or dark from
the input readout.  The overscan is taken from the \code{bias} member
of the input \code{pmReadout}, and can be turned off by supplying a
\code{NULL} value for \code{overscanOpts}.  A bias image is subtracted
if \code{bias} is non-\code{NULL}, and a dark frame is scaled (by
\code{CELL.DARKTIME}) and subtracted if \code{dark} is
non-\code{NULL}.
\begin{prototype}
bool pmSubtractBias(pmReadout *in,      ///< Input readout, to be overscan/bias/dark corrected
                    pmOverscanOptions *overscanOpts, ///< Options for overscan subtraction, or NULL
                    const pmReadout *bias, ///< Bias image to subtract, or NULL
                    const pmReadout *dark ///< Dark image to scale and subtract, or NULL
                   );
\end{prototype}


\subsubsection{Overscan subtraction}
\label{sec:overscan}

The options for performing the overscan subtraction are bundled in a
\code{pmOverscanOptions}:

\begin{datatype}
typedef struct {
    // Inputs
    bool single;                ///< Reduce all overscan regions to a single value?
    pmFit fitType;              ///< Type of fit to overscan
    unsigned int order;         ///< Order of polynomial, or number of spline pieces
    psStats *stat;              ///< Statistic to use when reducing the minor direction
    // Outputs
    psPolynomial1D *poly;       ///< Result of polynomial fit
    psSpline1D *spline;         ///< Result of spline fit
} pmOverscanOptions;
\end{datatype}

The overscan subtraction may be performed by reducing all overscan
regions to a single value (e.g., if there is no structure); or the
overscan may be fit perpendicular to the read direction (usually the
columns) with a particular functional form; or a single value may be
subtracted for each read/scan without fitting (if the structure defies
characterisation).  In any case, statistics are required to reduce
multiple values to a single value (either for the scan, or for the
entire overscan regions).

The \code{fitType} is an enumerated type which specifies the type of
fit to employed on the overscan vector:
\begin{datatype}
typedef enum {
    PM_FIT_NONE,                        ///< No fit
    PM_FIT_POLY_ORD,                    ///< Fit ordinary polynomial
    PM_FIT_POLY_CHEBY,                  ///< Fit Chebyshev polynomial
    PM_FIT_SPLINE                       ///< Fit cubic splines
} pmFit;
\end{datatype}

If \code{fitType} is \code{PM_FIT_NONE}, then the overscan vector is
subtracted from the image without fitting.  Otherwise, the overscan
vector is fit using the specified functional form, the fit is
subtracted from the image, and the \code{poly} or \code{spline} is
allocated and updated with the results of the fit.

The allocator for a \code{pmOverscanOptions} is:
\begin{prototype}
pmOverscanOptions *pmOverscanOptionsAlloc(bool single, ///< Reduce all overscan regions to a single value?
        pmFit fitType, ///< Type of fit to overscan
        unsigned int order, ///< Order of polynomial, or number of splines
        psStats *stat ///< Statistic to use
                                         );
\end{prototype}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Non-linearity}

We here specify two functions to perform the non-linearity correction,
since either (or both) might be used to specify the correction.

The first, \code{pmNonLinearityPolynomial}, applies a polynomial to
the flux of each pixel in the input image to determine the corrected
flux.  The second function, \code{pmNonLinearityLookup} performs
linear interpolation on the table (specified by two vectors) to
determine the corrected flux.

\begin{prototype}
pmReadout *pmNonLinearityPolynomial(pmReadout *in, ///< Input image, to correct
                                    const psPolynomial1D *coeff ///< Polynomial for non-linearity correction
                                   );
pmReadout *pmNonLinearityLookup(pmReadout *in, ///< Input image, to correct
                                const psVector *inFlux, ///< Table column with input fluxes
                                const psVector *outFlux ///< Table column with output fluxes
                               );
\end{prototype}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Shutter correction}

A mechanical shutter may not yield uniform exposure times as a
function of position on the detector.  The typical error consists of a
constant exposure-time offset relative to the requested value,
i.e. the exposure time is $T_o + \delta T(x,y)$.  The exposure error,
$\delta T$, may be measured with the following scheme.  Obtain a set
of exposures with different exposures times taken of the same
flat-field source; the source must be spatially stable between the
exposures, but need not have a stable amplitude.  For an illuminating
flux of intensity $F(x,y) = F_o f(x,y)$, the signal recorded by any
pixel in the detector is given by: $S(t,x,y) = F_o(t) f(x,y) (T_o +
\delta T(x,y))$ where $F_o$ is the $F_o(t)$ is the (variable) overall
intensity of the illuminating source and $f(x,y)$ is the spatial
illumination patter times the flat-field response.  Choose a reference
location in the image (e.g., the detector center) and divide by the
value of that region (i.e., mean or median):

\begin{eqnarray}
s(t,x,y) & = & S(t,x,y) / S(t,0,0) \\
s(t,x,y) & = & F_o(t) f(x,y) (T_o + \delta T(x,y)) / F_o(t) f(0,0) (T_o + \delta T(0,0)) \\
s(t,x,y) & = & f(x,y) (T_o + \delta T(x,y)) / f(0,0) (T_o + \delta T(0,0)) \\
\end{eqnarray}

We can absorb the term $f(0,0)$ into $f(x,y)$, as we have no
motivation for the scale of $f(x,y)$.  For any single pixel, over the
set of exposures, we thus need to solve for $\delta T(x,y)$, $\delta
T(0,0)$, and $f'(x,y)$ in the equation:

\begin{equation}
s(t,x,y) = f'(x,y) (T_o + dT(x,y)) / (T_o + \delta T(0,0))
\end{equation}

We avoid directly fitting these values, as the process would be a
non-linear least-squares problem for every pixel in the image, and
thus very time consuming.  There are linear options which may be used
instead.  First, as $T_o$ goes to a large value, $s$ approaches the
value of $f'(x,y)$.  Next, as $T_o$ goes to a very small value, $s$
approaches the value of $f'(x,y) \delta T(x,y)/ \delta T(0,0)$.
Finally, when s has the value of $f'(x,y)*(1 + \delta T(x,y)/ \delta
T(0,0))/2$, then $T_o$ has the value of $\delta T(0,0)$.  With data
points covering a reasonable dynamic range, we can solve for these
three values by interpolation and/or extrapolation.

To take the strategy one step further, we could use the above recipe
to obtain a guess for the three parameters and then apply non-linear
fitting to solve more accurately for the parameters.  If we limit this
operation to a handful of positions in the image (user defined, but
the obvious choice would be positions near the center, edges, and
corners), then we may determine a good value for $\delta T(0,0)$.
Since there is only one $\delta T(0,0)$ for the image, we can apply
the resulting measurement to the rest of the pixels in the image.  If
$\delta T(0,0)$ is not a free parameter, then the fitting process is
linear in terms of $\delta T(x,y)$ and $f'(x,y)$.

\subsubsection{Components to creating a shutter correction}

We define a shutter correction, applicable for a single pixel (or region):
\begin{datatype}
typedef struct {
    double scale;                       ///< The normalisation for an exposure, A(k)
    double offset;                      ///< The time offset, dTk
    double offref;                      ///< The reference time offset, dTo
} pmShutterCorrection;
\end{datatype}

The allocator is:
\begin{prototype}
pmShutterCorrection *pmShutterCorrectionAlloc();
\end{prototype}


The following function provides a guess for the shutter correction,
based on plot of counts vs exposure time.  It should be used before
doing the full non-linear fit, to get parameters close to the true.
The function assumes the \code{exptime} vector is sorted (ascending
order; longest is last) prior to input.
\begin{prototype}
pmShutterCorrection *pmShutterCorrectionGuess(const psVector *exptime, ///< Exposure times for each exposure
        const psVector *counts ///< Counts for each exposure
                                             );
\end{prototype}

A full non-linear fit may then be performed to the counts as a
function of exposure time.  The main purpose is to solve for the
reference time offset, so that future fits may be performed using
linear fitting with the reference time offset fixed.
\begin{prototype}
pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, ///< Exposure times for each exposure
        const psVector *counts, ///< Counts for each exposure
        const psVector *cntError, ///< Error in the counts
        const pmShutterCorrection *guess ///< Initial guess
                                               );
\end{prototype}

With the reference time offset fixed to a known value, a (much faster)
linear fit can be performed at each pixel to generate the shutter
correction.  This function performs the fit with iterative clipping,
if \code{nIter} $> 1$.
\begin{prototype}
pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, ///< Exposure times for each exposure
        const psVector *counts, ///< Counts for each exposure
        const psVector *cntError, ///< Error in the counts
        const psVector *mask, ///< Mask for each exposure
        float offref, ///< Reference time offset
        int nIter, ///< Number of iterations
        float rej, ///< Rejection threshold (sigma)
        psMaskType maskVal ///< Mask value
                                              );
\end{prototype}

\subsubsection{Creating a shutter correction}

The above steps are all performed using the following function.  Given
an array of images with known exposure times, this function measures
the shutter correction (our principal concern is for the time offset,
rather than the normalisation) by measuring the reference time offset
using the full non-linear fit for a small number of representative
regions (middle and corners), and then using that to perform a linear
fit to each pixel.
\begin{prototype}
psImage *pmShutterCorrectionMeasure(const psVector *exptimes, ///< Exposure times
                                    const psArray *images, ///< Input images
                                    const psArray *weights, ///< Weight images
                                    const psArray *masks, ///< Mask images
                                    unsigned int size, ///< Size of samples for statistics for non-linear fit
                                    psStatsOptions meanStat, ///< Statistic to use for mean
                                    psStatsOptions stdevStat, ///< Statistic to use for stdev
                                    int nIter, ///< Number of iterations
                                    float rej, ///< Rejection threshold (sigma)
                                    psMaskType maskVal ///< Mask value
                                   );
\end{prototype}

\subsubsection{Applying a shutter correction}

Given a shutter correction (a measurement of $\delta T$ for each
pixel), the following function applies this correction to an input
image.
\begin{prototype}
bool pmShutterCorrectionApply(pmReadout *readout, ///< Readout to which to apply shutter correction
                              const pmReadout *shutter ///< Shutter correction readout, with dT for each pixel
                             );
\end{prototype}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Flat-fielding}


The following function applies the flat-field calibration to the input
readout.  The relative offsets between the input and flat images is
determined from the readout row0,col0 and the CELL.X0 and CELL.Y0
concepts.  Normalisation of the flat is left as the responsibility of
the caller.  Non-positive pixels in the flat are masked, if there is a
mask present in the input readout.
\begin{prototype}
bool pmFlatField(pmReadout *in,         ///< Readout with input image
                 const pmReadout *flat  ///< Readout with flat image
                );
\end{prototype}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Masking}

The following function applies a mask readout to an input (science) readout.

Pixels marked as bad (according to the supplied mask value) within the
mask are marked as bad within the input image's mask.  If
\code{maskVal} is non-zero, all pixels in the mask have any of the
same bits sets as \code{maskVal} shall have the corresponding bits
raised.  If \code{maskVal} is zero, any zero pixels in the mask are
OR-ed with \code{PM_MASK_BAD}.  Position offsets (such as due to
trimming) between the input and mask are applied so that the same
pixels are referred to.  The science readout must already have a
supplied mask element (use e.g. \code{pmReadoutSetMask}).
\begin{prototype}
bool pmMaskBadPixels(pmReadout *input,  ///< Input science image
                     const pmReadout *mask, ///< Mask image to apply
                     psMaskType maskVal ///< Mask value to apply
                    );
\end{prototype}

\tbd{Masks for the charge traps need to be grown by the extent of the
OT convolution kernel.  For other pixel types, orthogonal transfer of
the flux in this pixel will not (necessarily) affect the flux in
neighbouring pixels.}


\subsection{Fringes}

Some images contain a signal caused by thin-film interference in the
device due to strong emission lines.  The resulting instrumental
effect consists of a pattern (the ``fringe pattern'') of bright and
dark bands corresponding to the constructive and destructive
interference of the emission lines.  In the case that a single
emission line causes the line structure, the resulting pattern can be
described by two independent parameters: First, the amplitude of the
emission line determines the overall amplitude of the pattern.
Second, the three-dimensional surface structure of the device
determines the shape of the pattern.  In a typical situation, the
device is illuminated by multiple emission lines, as well as a
continuum spectral source, which contributes to the overall light
detected by the device without following the fringe pattern.  The
relative intensities of the continuum background and the fringe
pattern depend on the device structure (thickness) and on the ratio of
the continuum and line emission fluxes. 

A simple approach to the fringe pattern is to subtract a master fringe
frame scaled by the amplitude of the fringe pattern.  The amplitude of
the fringe pattern is used both in the process of constructing the
master image and in scaling the master image when it is applied to
science image.  This is the method currently in use at CFHT and it
usually performs well.  However, the fringe signal can vary as the
emission lines in the atmosphere change, and the above method breaks
down unless different fringe images corresponding to different
atmospheric conditions are constructed.  

An alternative technique is to use multiple master fringe images at
the same time, each representing different atmospheric conditions.
The observed fringe frame can be considered as a linear combination of
different fringe patterns, depending on the relative strengths of the
lines active in creating each of the fringe masters.  It is not
critical that the fringe master images represent completely orthogonal
fringe patterns, they need only sample sufficiently different
conditions to provide a handle on the underlying fringe signals.  

We define a method of measuring the fringe pattern which is robust in
the presence of stars and which is fast.  We implement a varient on
the method used at CFHT in which the fringe pattern is mapped by a
series of points distributed across the image.  At CFHT, manual effort
is used to carefully define point pairs which correspond to peaks and
valleys of the fringe pattern.  We implement a different approach in
which the fringe points are randomly chosen across the image.  At each
point in the image, the median flux is measured in a box of specified
size.  A low-frequency spatial filter is then applied to these
measurements.  The resulting array of points and fluxes then
represents the strength of the fringe pattern on that image.  The
comparison between any two fringe images is then just a linear fit
between these fringe statistics vectors, as follows:
\[
S_i = C_0 + C_1 F_i
\]
where $S_i$ is the fringe statistic on the science image and $F_i$ is
the fringe statistic on the reference fringe image.  Extending this
logic to any number of reference fringe images results in the
following relationship:
\[
S_i = C_0 + \sum_j C_j F_j
\]

In order to correct a single science image, the collection of fringe
statistics ($S_i$) are used to measure the coefficients $C_0$, $C_j$.
The linear combination of the reference fringe images is then used to
build a master image which is subtracted from the science image.  The
following structures and functions implement the above concepts.


\subsubsection{Fringe regions}

Fringes are measured within a box of size \code{dX,dY}.  A large scale
smoothing is performed by subtracting the background within large
divisions of the image.  The coordinates of the fringe points and the
mask are originally set to \code{NULL}, which means that they will be
generated when required.

\begin{datatype}
typedef struct {
    psU32 nRequested;                   // Number of fringe points selected
    psU32 nAccepted;                    // Number of fringe points not masked
    psU32 dX;                           // Median box half-width
    psU32 dY;                           // Median box half-height
    psU32 nX;                           // Number of large-scale smoothing divisions in x (col)
    psU32 nY;                           // Number of large-scale smoothing divisions in y (row)
    psVector *x;                        // Fringe point coordinates (col), or NULL
    psVector *y;                        // Fringe point coordinates (row), or NULL
    psVector *mask;                     // Fringe point on/off mask, or NULL
} pmFringeRegions;
\end{datatype}

The allocator is:
\begin{prototype}
pmFringeRegions *pmFringeRegionsAlloc (int nPts, ///< Number of fringe points to create
                                       int dX, ///< Half-width of fringe boxes
                                       int dY, ///< Half-height of fringe boxes
                                       int nX, ///< Smoothing scale in x
                                       int nY ///< Smoothing scale in y
                                      );
\end{prototype}

The following function generates fringe points randomly over an image.
No effort is made to avoid masked regions (indeed, the function knows
nothing about masks).  If the random number generator is \code{NULL},
then a new one will be used.
\begin{prototype}
bool pmFringeRegionsCreatePoints(pmFringeRegions *fringe, ///< Fringe regions to generate
                                 const psImage *image, ///< Image for the regions (defines the size)
                                 psRandom *random ///< Random number generator, or NULL
                                );
\end{prototype}

Fringe regions may be written to and read from FITS files with the following functions.

\code{pmFringeRegionsWriteFits} writes fringe regions to a FITS file,
with the given extension name.  The header is supplemented with scalar
values \code{dX}, \code{dY}, \code{nX} and \code{nY} (as
\code{PSFRNGDX}, \code{PSFRNGDY}, \code{PSFRNGNX}, \code{PSFRNGNY})
from the fringe regions, while the fringe coordinates and mask are
written as a FITS table (as \code{x}, \code{y}, \code{mask}).

\code{pmFringeRegionsReadFits} reads fringe regions from a FITS file,
at the given extension name.  The scalars are retrieved from the
header, while the table provides the fringe coordinates and mask.

\begin{prototype}
bool pmFringeRegionsWriteFits(psFits *fits, ///< Output FITS file
                              psMetadata *header, ///< Additional headers to write, or NULL
                              const pmFringeRegions *regions, ///< Regions to write
                              const char *extname ///< Extension name, or NULL
                             );
pmFringeRegions *pmFringeRegionsReadFits(psMetadata *header, ///< Header to read, or NULL
        const psFits *fits, ///< Input FITS file
        const char *extname ///< Extension name, or NULL
                                        );
\end{prototype}

\subsubsection{Fringe Statistics}

Fringe measurements for a particular image are contained in the
following structure, which contains the median and standard deviation
for each of the fringe regions.
\begin{datatype}
typedef struct {
    pmFringeRegions *regions;           ///< Fringe regions
    psVector *f;                        ///< Fringe point median
    psVector *df;                       ///< Fringe point stdev
} pmFringeStats;
\end{datatype}

The allocator is:
\begin{prototype}
pmFringeStats *pmFringeStatsAlloc(pmFringeRegions *regions // The fringe regions which will be measured
                                 );
\end{prototype}

The following function measures the fringe statistics (median and
standard deviation) at each of the fringe regions for an image
(supplied within a readout).  If the fringe regions are undefined,
they are generated.
\begin{prototype}
pmFringeStats *pmFringeStatsMeasure(pmFringeRegions *fringe, ///< Fringe regions at which to measure
                                    const pmReadout *readout, ///< Readout for which to measure
                                    psMaskType maskVal ///< Mask value for image
                                   );
\end{prototype}

Like the fringe regions, the fringe statistics may be written to and
read from a FITS file, allowing rapid re-generation of the fringe
measurements for the calibration files.

\code{pmFringeStatsWriteFits} writes the fringe statistics for an
image to a FITS table, with the given extension name.  The median and
standard deviation measurements are written as a FITS table (as
\code{f} and \code{df}).

\code{pmFringeStatsReadFits} reads the fringe stats from a FITS table,
at the given extension name.  It is assumed that the fringe
measurements correspond to the regions provided.

\begin{prototype}
bool pmFringeStatsWriteFits(psFits *fits, ///< FITS file to which to write
                            psMetadata *header, ///< Additional headers to write, or NULL
                            const pmFringeStats *fringe, ///< Fringe statistics to be written
                            const char *extname ///< Extension name for table
                           );
pmFringeStats *pmFringeStatsReadFits(psMetadata *header, ///< Header to read, or NULL
                                     const psFits *fits, ///< FITS file from which to read
                                     const char *extname, ///< Extension name to read
                                     pmFringeRegions *regions ///< Corresponding regions
                                    );
\end{prototype}

Since each readout of each chip must be measured separately (so as to
avoid any gaps between the cells, as in the case for GPC), while the
fit must be performed with all the readouts belonging to a chip (in
order to get a secure measurement of the fringe amplitudes), we
provide the following function, which concatenates the fringe stats
for several readouts into a single fringe stats.  The function
generates a new \code{pmFringeStats} from concatenating those in the
array.  The corresponding \code{pmFringeRegions} is also generated.
\begin{prototype}
pmFringeStats *pmFringeStatsConcatenate(const psArray *fringes, ///< Array of pmFringeStats for the readouts
                                        const psVector *x0, ///< Offset in x for the readout
                                        const psVector *y0 ///< Offset in y for the readout
                                       );
\end{prototype}


\subsubsection{Fringe scales}

The following structure contains the fringe solution --- the
coefficients of the linear system of library fringe images that best
fit the science image.

\begin{datatype}
typedef struct {
    int nFringeFrames;                  ///< Number of fringe frames
    psVector *coeff;                    ///< Fringe coefficients; size = nFringeFrames
    psVector *coeffErr;                 ///< Error in fringe coefficients; size = nFringeFrames
} pmFringeScale;
\end{datatype}

The corresponding allocator is:
\begin{prototype}
pmFringeScale *pmFringeScaleAlloc(int nFringeFrames ///< Number of fringe frames
                                 );
\end{prototype}

The following function measures the fringe scales for the library
fringe images, solving the linear system that best fits the fringe
measurements made on the science image.  Given a fringe measurement
for a science image, and an array of template fringe measurements, the
function measures the contribution of each of the templates to the
input.  Rejection is performed on the fringe regions, to weed out
stars etc.
\begin{prototype}
pmFringeScale *pmFringeScaleMeasure(pmFringeStats *science, ///< Fringe measurements from science image
                                    psArray *fringes, ///< Array of fringe measurements from templates
                                    float rej, ///< Rejection threshold (in standard deviations)
                                    unsigned int nIter, ///< Maximum number of iterations
                                    float keepFrac ///< Minimum fraction of regions to keep
                                   );
\end{prototype}


\subsubsection{Fringe correction}

The following function solves for and applies the fringe correction.
It is a wrapper around each of the fringe correction components to
measure the fringe points, solve for the fringe correction, and apply
the fringe correction.  The input fringe images are modified (scaled
by the solution coefficients in order to correct the science image).
The function returns the summed fringe image, and corrects the input
science image.

\begin{prototype}
psImage *pmFringeCorrect(pmReadout *in, ///< Input science image
                         pmFringeRegions *fringes, ///< The fringe regions used
                         psArray *fringeImages, ///< Fringe template images to use in correction
                         psArray *fringeStats, ///< Fringe stats (for templates) to use in correction
                         psMaskType maskVal, ///< Value to mask for science image
                         float rej,     ///< Rejection threshold, for pmFringeScaleMeasure
                         unsigned int nIter, ///< Maximum number of iterations, for pmFringeScaleMeasure
                         float keepFrac ///< Minimum fraction of regions to keep, for pmFringeScaleMeasure
                        );
\end{prototype}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Subtract sky}

\tbd{This has been deferred.}

Given an input image, a polynomial or spline specifying the order of a
desired fit, a binning factor and statistics to use for the binning,
along with a clipping level, \code{pmSubtractSky} shall fit and
subtract a model for the background of the image.  The API shall be
the following:
\begin{prototype}
pmReadout *pmSubtractSky(pmReadout *in, psPolynomial2D *poly, psImage *mask, psU8 maskVal, 
                         int binFactor, psStats *stats, float clipSD);
\end{prototype}

Note that the input image, \code{in}, shall be subtracted in-place.
The function shall return the subtracted image, and also update the
polynomial, Chebyshev or spline specified by \code{fitSpec}, to hold
the coefficients used in the subtraction.

The polynomial, \code{poly}, specifies the order of the polynomial,
and on return shall contain the coefficients of the fit.  If
\code{poly} is \code{NULL}, then no fit shall be performed, and the
function shall generate a warning and return.

When fitting the polynomial, the function shall first bin the input
image by \code{binFactor} in order to reduce the required processing
time.  In the binning, pixels in the \code{mask} (if non-\code{NULL})
which satisfy the \code{maskVal} shall be excluded.  The statistic to
use in this binning is specified by \code{stat}.  \code{stat} is of
type \code{psStats} instead of simply \code{psStatsOptions} so that
clipping levels may be specified, if desired.  In the event that
multiple options are specified by \code{stats}, a warning shall be
generated, and the option with the highest priority shall be used,
according to the following priority order: \code{PS_STAT_SAMPLE_MEAN},
\code{PS_STAT_SAMPLE_MEDIAN}, \code{PS_STAT_CLIPPED_MEAN},
\code{PS_STAT_ROBUST_MEAN}, \code{PS_STAT_ROBUST_MEDIAN},
\code{PS_STAT_ROBUST_MODE}.  If the \code{binFactor} is non-positive,
or \code{stats} is \code{NULL} or fails to specify an option, a
warning shall be generated, and the fit shall be performed on the
entire image.

Binned pixels deviating more than \code{clipSD} standard deviations
from the mean of the binned pixels shall be clipped in a single
clipping iteration before polynomial fitting.  These pixels may be
interpolated over, or may be simply ignored in the fitting, according
to the choice of algorithm.  If the \code{clipSD} is non-positive,
then the function shall generate a warning and not perform any
clipping.

The \code{mask} shall be of type U8, and the input image,
\code{in}, must be of type F32.

This operation acts only on the region of the readout specified by
\code{CELL.TRIMSEC}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Paper Trail}

\tbd{This is not implemented yet.}

The elements of the focal plane hierarchy each contain an
\code{analysis} member, intended to log the results of the detrend
tasks.  The detrend tasks shall add to the \code{analysis} members as
follows:

\begin{itemize}
\item \code{pmMaskBadPixels}:
  \begin{itemize}
  \item \code{MASK.DONE} (STR): The time at which masking was
    completed.
  \item \code{MASK.SAT} (S32): The number of saturated pixels masked
    in the image
  \item \code{MASK.SAT.GROW} (S32): The number of additional pixels
    masked by growing the saturated pixels.
  \item \code{MASK.BAD} (S32): The number of pixels masked in the
    image
  \item \code{MASK.BAD.GROW} (S32): The number of additional pixels
    masked by growing the specified bad pixels.
  \end{itemize}
\item \code{pmNonLinearityPolynomial} and \code{pmNonLinearityLookup}:
  \begin{itemize}
  \item \code{NONLIN.DONE} (STR): The time at which the non-linearity
    correction was completed.
  \item \code{NONLIN.POLY} (STR): The polynomial coefficients used (if
    applicable).
  \item \code{NONLIN.LOOKUP} (STR): The filename for the lookup table
    (if applicable).
  \end{itemize}
\item \code{pmSubtractBias}:
  \begin{itemize}
  \item \code{BIAS.DONE} (STR): The time at which the bias-subtraction
    was completed.
  \item \code{BIAS.OVERSCAN.AXIS} (STR): Overscan axis used.
  \item \code{BIAS.OVERSCAN.FIT.TYPE} (STR): Fit type applied to
    overscan.
  \item \code{BIAS.OVERSCAN.FIT.COEFF} (STR): Coefficients of overscan
    fit.
  \item \code{BIAS.OVERSCAN.REGION} (STR): Overscan regions (from
    \code{x0,y0,numCols,numRows}).
  \item \code{BIAS.OVERSCAN.BIN} (S32): Number of pixels per bin used
    in overscan.
  \item \code{BIAS.OVERSCAN.MEAN} (F32): The mean of the binned
    overscan pixels after subtracting the fit.
  \item \code{BIAS.OVERSCAN.SD} (F32): The standard deviation of the
    binned overscan pixels after subtracting the fit.
  \end{itemize}
\item \code{pmFlatField}:
  \begin{itemize}
  \item \code{FLAT.DONE} (STR): The time at which the flat-fielding
    was completed.
  \item \code{FLAT.BAD} (S32): Number of non-positive flat-field
    pixels.
  \end{itemize}
\end{itemize}

To be added by higher-levels:
\begin{itemize}
\item \code{BIAS.NAME} (STR): Name of bias image
\item \code{DARK.NAME} (STR): Name of dark image
\item \code{FLAT.NAME} (STR): Name of flat image
\item \code{MASK.NAME} (STR): Name of mask image
\end{itemize}

\subsection{Detrend Lookups}

When it comes time to perform a detrend operation on an image, it is
necessary to determine {\em which} detrend image should be used.  The
Pan-STARRS Image Processing Pipeline uses the concept of a detrend
image database table, or set of tables (part of the Metadata
Database), to store the known master detrend images.  These tables can
be accessed though the basic query functions specified for the master
detrend database.  To simplify the interaction for the case of the
detrend images, the following function allows the user to explicitly
search the detrend database table or tables for detrend images which
satisfy a set of characteristics.

\begin{prototype}
psArray *pmDetrendLookup (psMetadata *constraints, psMetadata *tableDefs);
\end{prototype}
This function accepts a metadata structure which restricts the
selected detrend images.  This metadata structure may contain any of
the following entries:
\begin{verbatim}
TYPE        type of detrend data (eg, flat, bias) 
CAMERA      name of desired camera (eg, GPC, MEGACAM)
CHIP        chip identifier (eg., ccd00)
FILTERNAME  name of specific filter hardware (eg, r.GPC01)
FILTERTYPE  conceptual name of filter (eg., r)
TIME_MIN    lower bound on valid time range 
TIME_MAX    upper bound on valid time range 
LABEL       match the entry label
RECIPE      recipe used to build detrend image
EXPTIME     exposure time 
AIRMASS     airmass 
\end{verbatim}
Any detrend images which match the provided constraints are returned
as an array of \code{psMetadata} elements corresponding to the columns
of the detrend database table.  The additional input parameter
specifies additional information to define the detrend database
tables.  This may include the access information (IP, Username,
Password), as well as names for the table and the columns which
correspond to the constraint names.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Detrend Creation}

In the detrend creation process, a collection of raw images are
combined to produce a clean, high-quality master image for correcting
the effect of interest.  The input images may potentially be processed
and scaled in some way.  The resulting output images may be to be
re-scaled to have a consistent signal for all chips in the mosaic.
The simplest example is the construction of a bias image, in the case
where there is signficant 2-D bias structure.  In this case, the input
raw bias images are probably combined without any additional
processing.  In another example, flat-field image must be
bias-corrected and scaled to a consistent normalization before being
combined, and the flat-field images from the different chips must be
normalized so that each chip will be flattened consistently across the
mosaic.  A complex example is the fringe pattern, in which the input
images must be bias-corrected and flattened, and the resulting images
must be scaled by the amplitude of the fringe pattern on each image,
rather than by the average flux level.  In this section, we define the
tools necessary to perform the detrend creation process.

\subsection{Image Stacking}

A basic operation in generating the master detrend images is using a
stack of many input images of a particular type and combining them,
with perhaps some additional scaling, in order to build up
signal-to-noise and to reject deviant pixel.  For this, we require a
general purpose image combination module.  We forsee this module as
only acting upon data from the same detector, and so each input image
will have the same noise characteristics.

The following structure specifies how the combination of a stack of
pixels is performed.  The particular statistic specified by
\code{combine} is used to combine each stack of pixels from the input
images.  If the \code{combine} specifies clipped fitting, the
\code{iter} and \code{rej} members provide the number of iterations
and the rejection threshold for clipping.

The \code{maskVal} specifies mask values for pixels to be
masked out of the combination.

After masking, but before performing the combination, the highest
\code{fracHigh} fraction and lowest \code{fracLow} fraction of pixels
in the stack are immediately rejected, unless this would leave less
than \code{nKeep} pixels in the stack, in which case no immediate
rejection is performed.

\begin{datatype}
typedef struct {
    psStatsOptions combine;             ///< Statistic to use when performing the combination
    psMaskType maskVal;                 ///< Mask value
    int nKeep;                          ///< Mimimum number of pixels to keep
    float fracHigh;                     ///< Fraction of high pixels to immediately throw
    float fracLow;                      ///< Fraction of low pixels to immediately throw
    int iter;                           ///< Number of iterations for clipping (for CLIPPED_MEAN only)
    float rej;                          ///< Rejection threshould for clipping (for CLIPPED_MEAN only)
} pmCombineParams;
\end{datatype}

The allocator is:
\begin{prototype}
pmCombineParams *pmCombineParamsAlloc(psStatsOptions statsOptions ///< Statistic to use for combination
                                     );
\end{prototype}

The combination is performed by the following function:
\begin{prototype}
bool pmReadoutCombine(pmReadout *output,///< Output readout; altered and returned
                      const psArray *inputs,  ///< Array of input readouts (F32 image and weight, U8 mask)
                      const psVector *zero, ///< Zero corrections to subtract from input, or NULL
                      const psVector *scale, ///< Scale corrections to divide into input, or NULL
                      const pmCombineParams *params ///< Combination parameters
                     );
\end{prototype}

\code{pmReadoutCombine} combines input images pixel by pixel --- for
each pixel of the output image, a stack of contributing input pixels
is formed and combined, using the mask and weights if supplied in the
input readouts.

The \code{output} readout is extended in size, if necessary, to
contain the input pixels.

The \code{inputs} array contains the \code{pmReadout}s to be combined.
The images contained within the \code{pmReadout}s need not all be of
the same size, but the module will take into account the offsets from
the corner of the detector when comparing pixels, so that it is the
same \textit{physical} pixels that are combined.

If the \code{zero} vector is non-\code{NULL} then the appropriate
values are subtracted from the \code{inputs} before rejection is
performed.  If the \code{scale} vector is non-\code{NULL}, then the
appropriate values multiply the \code{inputs} before rejection is
performed.


\subsection{Flat-field Re-Normalization}

Consider a collection of $N_i$ flat-field images obtained with a
mosaic camera consisting of $N_j$ chips.  Each image is exposed to an
illumination source which should be a uniform surface
brightness\footnote{This is likely a false assumption: the
illumination source likely has spatial variations.  However, for the
purposes of this discussion, it only matters that such spatial
variations scale consistently as a function of illumination intensity.
The spatial errors are corrected by the photometric flat-field
correction technique (eg., Magnier \& Cuillandre 2004).}  Two factors
determine the actual measured flux level (in Digital Numbers) on each
of the chips in each image: the gain of each chip ($\mbox{gain}_j$)
and the flux level from the illumination source ($\mbox{source}_i$).
When the images are combined, the input images must be scaled so that
the flux levels can be consistently compared.  After combining the
collection of images, it is necessary to determine an appropriate
re-normalization for the resulting flat-field images.  In effect, the
individual chips must be adjusted so that the master flat-field image
has a flux level which varies from chip to chip in proportion to the
actual chip gain.  In this case, if a uniform illumination source
illuminates the mosaic, the resulting flux levels will be corrected by
the flat-field to a single, consistent flux level. 

In order to determine the correct relative scaling between the
devices, it is thus necessary to know the individual chip gains, or at
least the gain ratios.  A typical technique scaled all chips relative
to a reference chip, or by a statistic measured for the complete
collection.  These techniques fail if the input collection of images
does not always consist of the same set of chips; for the GPC on
Pan-STARRS, we must expect that individual cells or even chips may be
disabled on a frequent basis, so our algorithms must not be limited by
the assumption that all chips are available in all images.  We
therefore define the following algorithm to measure the relative chip
gains for a collection of input flat-field images, each with a
measured flux $\mbox{flux}_{i,j}$.  We want to solve for the chip
gains and the source illumination fluxes which would make the best
prediction of the measured input image fluxes:
\[
\mbox{flux}^{\rm pred}_{i,j} = \mbox{gain}_j \times \mbox{source}_i
\]
This relationship is easiest to determine if we take the logarithm of
both sides of the equation:
\[
M^{\rm pred}_{i,j} = G_j + S_i
\]
where $M^{\rm pred}_{i,j} = \log (\mbox{flux}^{\rm pred}_{i,j})$, $G_j
= \log (\mbox{gain}_j)$, and $S_i = \log (\mbox{source}_i)$.  We can
then write the chi-square which we want to minimize as:
\[
\chi^2 = \sum_{i,j} (M^{\rm obs}_{i,j} - G_j - S_i)^2
\]
where we ignore the weights of the different measured flux levels.
Taking the derivatives with respect to the parameters of interest
($G_j, S_i$), and setting them to 0, we determine the following set of
equations which must be solved:
\[
G_j \times N_i = \sum_i M^{\rm obs}_{i,j} - \sum_i S_i
\]
\[
S_j \times N_j = \sum_j M^{\rm obs}_{i,j} - \sum_j G_j \\
\]
This set of equations can be solved iteratively, starting from the
assumption that all chip gains are 1.0, ($G_j = 0$), or by supplying
a guess for the chip gains.  The result of this analysis is the
measured chip gains and the measured source illumination levels for
each of the input flat-field images.  The chip gains can then be used
to modify the flux levels on the master flat-field images.

The following function solves this system.  The matrix of background
measurements ($M_{i,j}$) contains the background for the flat fields
used in the combination, as a function of exposure (rows) and chip
(columns).  The exposure fluxes and chip gains are modified upon
return with the solved values.  The function returns true if the
solution converged.  An initial guess for the chip gains might be
helpful, but is not necessary.
\begin{prototype}
bool pmFlatNormalize(psVector **expFluxesPtr, ///< Flux in each exposure, or NULL; modified
                     psVector **chipGainsPtr, ///< Initial guess of the chip gains or NULL; modified
                     const psImage *bgMatrix ///< Background measurements: rows are exposures, cols are chips
                    );
\end{prototype}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Objects on Images}

\tbd{EAM to update, if necessary.}

\subsection{Overview}

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

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

In the discussion that follows, images are of type F32 and masks are
of type U8.

\subsection{Structures to Describe Sources}

In the object analysis process, we will use specific mask values to
mark the image pixels.  The following structure defines the relevant
mask values.
\begin{datatype}
typedef enum {
    PSPHOT_MASK_CLEAR     = 0x00,
    PSPHOT_MASK_INVALID   = 0x01,
    PSPHOT_MASK_SATURATED = 0x02,
    PSPHOT_MASK_MARKED    = 0x08,
} psphotMaskValues;
\end{datatype}

\subsubsection{pmSource and pmPeak}

We define the following structure to represent a single source
detected in a single image.  
\begin{datatype}
typedef struct {
  pmPeak *peak;            // description of peak pixel
  psImage *pixels;         // rectangular region including object pixels
  psImage *weight;         // Image variance.
  psImage *mask;           // Mask which marks pixels associated with objects.
  pmMoments *moments;      // basic moments measure for the object
  pmModel *modelPSF;       // PSF model parameters and type
  pmModel *modelEXT;       // FLT model parameters and type
  pmSourceType type;       // Best identification of object
  pmSourceMode mode;       // flags describing the model quality 
  psArray *blends;         // array of other sources blended with this source
  float apMag;             // measured aperture magnitude for source
  float fitMag;            // measured model magnitude for source
  psRegion region;         // area on image covered by selected pixels
} pmSource;
\end{datatype}

A source has the capacity for several types of measurements.  The
simplest measurement of a source is the location and flux of the peak
pixel associated with the source:
\begin{datatype}
typedef struct {
  int x;                   // x-coordinate of peak pixel
  int y;                   // y-coordinate of peak pixel
  float counts;            // value of peak pixel (above sky?)
  pmPeakType class;        // description of peak
} pmPeak;
\end{datatype}

A peak pixel may have several features which may be determined when
the peak is found or measured.  These are specified by the
\code{pmPeakType} enum.  \code{PM_PEAK_LONE} represents a single pixel
which is higher than its 8 immediate neighbors.  The
\code{PM_PEAK_EDGE} represents a peak pixel which touching the image
edge.  The \code{PM_PEAK_FLAT} represents a peak pixel which has more
than a specific number of neighbors at the same value, within some
tolerance:
\begin{datatype}
typedef enum {
  PM_PEAK_LONE,             // isolated peak
  PM_PEAK_EDGE,             // peak on edge
  PM_PEAK_FLAT              // peak has equal-value neighbors
  PM_PEAK_UNDEF             // Undefined.
} pmPeakType; 
\end{datatype}

\subsubsection{pmMoments and source description}

The pixels which contain the source are specified with the
\code{psImage *pixels} element, a subimage of the image being
analysed.  Similarly, the \code{mask} element is a subimage of the
corresponding mask image and the \code{weight} element is a subimage
of the corresponding weight image (image varience).  Since these are
subimages, a collection of many objects may include overlapping
pixels; care must be taken that pixel manipulations for one source do
not unintentionally interfere with the other source pixels.  The
\code{mask} may be used to exclude any pixels which are not considered
part of the source.  Along with these pixel structures, we include the
\code{psRegion region} element which defines the boundaries of the
current associated subimages. 

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{datatype}
typedef struct {
  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)
  float SN;                 // approx signal-to-noise
  int   nPixels;            // number of pixels used
} pmMoments;
\end{datatype}

A collection of object moment measurements can be used to determine
approximate object classes.  The key to this analysis is the location
and statistics (in the second-moment plane, $\sigma_x$ vs $\sigma_y$)
of the group of objects which are likely PSF objects.  We define the
following structure to identify the location and size of the psf clump
in the second-moment plane.
\begin{datatype}
typedef struct {
    float X;
    float dX;
    float Y;
    float dY;
} pmPSFClump;
\end{datatype}

A given source may be identified as most-likely to be one of several
source types.  The \code{pmSource} entry \code{pmSourceType} defines
the current best-guess for this source.  

\begin{datatype}
typedef enum {
    PM_SOURCE_UNKNOWN,                  ///< no guess yet made
    PM_SOURCE_DEFECT,                   ///< a cosmic-ray
    PM_SOURCE_SATURATED,                ///< random saturated pixels
    PM_SOURCE_STAR,                     ///< a good-quality star
    PM_SOURCE_EXTENDED,                 ///< an extended object (eg, galaxy)
} pmSourceType;
\end{datatype}

The related element, \code{pmSourceMode mode}, holds a collection of flags which
are used to indicate the status of the analysis for a source.  These
are defined below:
\begin{datatype} 
typedef enum {
    PM_SOURCE_DEFAULT    = 0x0000, ///< no flags are set
    PM_SOURCE_PSFMODEL   = 0x0001, ///< flags refer to the PSF model
    PM_SOURCE_EXTMODEL   = 0x0002, ///< flags refer to the EXT model
    PM_SOURCE_SUBTRACTED = 0x0004, ///< the model has been subtracted from the image
    PM_SOURCE_FITTED     = 0x0008, ///< the source has been fitted with a model
    PM_SOURCE_FAIL       = 0x0010, ///< the model fit failed
    PM_SOURCE_POOR       = 0x0020, ///< the model fit was poor (low S/N, etc)
    PM_SOURCE_PAIR       = 0x0040, ///< the model fit is one of a paired source
    PM_SOURCE_PSFSTAR    = 0x0080, ///< the source was used to construct the image PSF model
    PM_SOURCE_SATSTAR    = 0x0100, ///< the source is saturated
    PM_SOURCE_BLEND      = 0x0200, ///< the source is a blend with another source
    PM_SOURCE_LINEAR     = 0x0400, ///< the source was fitted with the linear PSF model
    PM_SOURCE_TEMPSUB    = 0x0800, ///< the source has been subtracted, but should be replaced
} pmSourceMode;
\end{datatype}

\subsubsection{pmModel Source Model and Abstraction} 

An object's flux distribution may be modeled 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
identified by code \code{type}, dynamically assigned based on the
available models (see below).  We discuss the details of these models
in section~\ref{ObjectModels}.  The model parameters have 4 special
elements.  The first four elements represent aspects of the source
which are not specified by the image PSF, even for point sources.  
These consist of, in order:
\begin{itemize}
\item the local sky
\item the object normalization
\item the x-coordinate
\item the y-coordinate
\end{itemize}

\tbd{should be include utility pointers to these parameters so that
  functions do not need to know the parameter sequence?}

The structure which carries the information about a given source model
is defined below:
\begin{datatype}
typedef struct {
  pmModelType type;         // model to be used
  psVector *params;         // parameter values
  psVector *dparams;        // parameter errors
  psF32 chisq;              // fit chisq
  psS32 nDOF;               // number of degrees of freedom
  psS32 nIter;              // number of iterations
  pmModelStatus status;     // fit status
  float radius;             // fit radius actually used
} pmModel;
\end{datatype}

The \code{status} element carries the resulting success/failure status
of an attempt to fit the model to the source:
\begin{datatype}
typedef enum {
    PM_MODEL_UNTRIED,               ///< model fit not yet attempted
    PM_MODEL_SUCCESS,               ///< model fit succeeded
    PM_MODEL_NONCONVERGE,           ///< model fit did not converge
    PM_MODEL_OFFIMAGE,              ///< model fit drove out of range
    PM_MODEL_BADARGS                ///< model fit called with invalid args
} pmModelStatus;
\end{datatype}

We distinguish several ways in which an analytical model may be
applied to a source.  The PSF model represents the best fit of the
image PSF to the specific object.  In this case, the PSF-dependent
parameters are specified for the object by the PSF, not by the fit.
The EXT model represents the best fit of the given model to the
object, with all parameters floating in the fit.  Such a model would
typically be used to represent and extended object, hence the
abbreviation EXT.  In some circumstances, a source may be fitted with
a PSF model in which the position is held fixed, and not allowed to
vary in the model fitting process.  We identify such a model as FIX.
Finally, we allow for the case in which two nearly-merged PSFs are
fitted with a single 2-PSF model.  We identify such a model as DBL.
The \code{pmSource} structure contains a pointer to both a PSF and an
EXT model, allowing any source to carry information about both
possible fitting modes \tbd{not clear that we actually use this
information; we might be better off simply distinguishing with one of
the pmSourceMode flags}.  The value of the model at a specific
coordinate can be determined by calling the function:
\begin{prototype}
psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row);
\end{prototype}
For this function, the values of \code{col,row} are in the
\code{image} coordinates, which may be a subimage, while the reference
coordinate for the model is in the parent image coordinates.

In the \code{pmSource} structure, the elements \code{apMag} and
\code{fitMag} are used to carry the measured magnitude of the source
determined either from aperture photometry or from the integral of the
fitted model function.  The element \code{blends} is used to carry
pointers to the collection of sources which were found to be blended
with this source.  Only the primary source of a blend group carries
this information.%%% (see Section~\ref{blends}).

Every model instance belongs to a class of models, defined by the
value of the \code{pmModelType type} entry.  Various functions need
access to information about each of the models.  Some of this
information varies from model to model, and may depend on the current
parameter values or other data quantities.  In order to keep the code
from requiring the information about each model to be coded into the
low-level fitting routines, we define a collection of functions which
allow us to abstract this type of model-dependent information.  These
generic functions take the model type and return the corresponding
function pointer for the specified model.  Each
model is defined by creating this collection of specific functions,
and placing them in a single file for each model.  We define the
following structure to carry the collection of information about the
models. 

\begin{datatype}
typedef struct {
    char *name;
    int nParams;
    pmModelFunc          modelFunc;
    pmModelFlux          modelFlux;
    pmModelRadius        modelRadius;
    pmModelLimits        modelLimits;
    pmModelGuessFunc     modelGuessFunc;
    pmModelFromPSFFunc   modelFromPSFFunc;
    pmModelFitStatusFunc modelFitStatusFunc;
} pmModelGroup;
\end{datatype}

Each entry in the \code{pmModelGroup} defines the information needed
by the system to specify a model.  The function types define above are
\begin{prototype}
typedef psMinimizeLMChi2Func pmModelFunc;
typedef psF64 (*pmModelFlux)(const psVector *params);
typedef psF64 (*pmModelRadius)(const psVector *params, double flux);
typedef bool (*pmModelLimits)(psVector **beta_lim, psVector **params_min, psVector **params_max);
typedef bool (*pmModelGuessFunc)(pmModel *model, pmSource *source);
typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf);
typedef bool (*pmModelFitStatusFunc)(pmModel *model);
\end{prototype}

Each of these functions is found for a given model by calling the
corresponding lookup function:
\begin{prototype}
pmModelFunc          pmModelFunc_GetFunction (pmModelType type);
pmModelFlux          pmModelFlux_GetFunction (pmModelType type);
pmModelRadius        pmModelRadius_GetFunction (pmModelType type);
pmModelLimits        pmModelLimits_GetFunction (pmModelType type);
pmModelGuessFunc     pmModelGuessFunc_GetFunction (pmModelType type);
pmModelFromPSFFunc   pmModelFromPSFFunc_GetFunction (pmModelType type);
pmModelFitStatusFunc pmModelFitStatusFunc_GetFunction (pmModelType type);
\end{prototype}

\code{pmModelFunc} is the function used to determine the value of the
model at a specific coordinate, and is the one used by
\code{psMinimizeLMChi2}.  

\code{pmModelFlux} returns the total integrated flux for the given
input parameters.

\code{pmModelRadius} returns the scaling radius at which the flux of
the model matches the specified flux.  This presumes that the model is
a function of an elliptical contour.  

\code{pmModelLimits} sets the parameter limit vectors for the
function.

\code{pmModelGuessFunc} generates an initial guess for the model based
on the provided source statistics (moments and pixel values as
needed).

\code{pmModelFromPSFFunc} takes as input a representation of the psf
and a value for the model and fills in the PSF parameters of the
model.  The input primarily relies upon the centroid coordinates of
the input model, thought the normalization may potentially be used.

\code{pmModelFitStatusFunc} returns a true or false values based on
the success or failure of a model fit.  the success is determined by
quantities such as the chisq or the signal-to-noise.

In addition, the following functions are useful for interacting with
the collection of models:
\begin{prototype}
int                  pmModelParameterCount (pmModelType type);
\end{prototype}
This function returns the number of parameters used by the listed
function.

\begin{prototype}
char                *pmModelGetType (pmModelType type);
pmModelType          pmModelSetType (char *name);
\end{prototype}
These two functions provide translations between the user-space model
names and the internal model type codes.  The model type codes are not
necessarily maintained between compilations of the program; the name
should be used to transfer models between programs or systems.

\subsubsection{pmGrowthCurve}

When the photometry of source is measured in a fixed aperture, there
is always a fraction of the source light which falls outside of the
aperture.  The resulting aperture magnitude is thus larger (ie,
fainter) than the actual source.  As the aperture is increased, the
amount of loss decreases and the measured magnitude increases.  This
trend is the curve of growth for the source.  We use the following
structure to carry information about the curve of growth.  We use the
PSF model to measure the curve of growth for an image.  

\begin{datatype}
typedef struct {
    psVector *radius;
    psVector *apMag;
    psF32 refRadius;
    psF32 maxRadius;
    psF32 fitMag;
    psF32 apRef;   // apMag[refRadius]
    psF32 apLoss;  // fitMag - apRef
} pmGrowthCurve;
\end{datatype}
In this structure, \code{radius} is a monotonically increasing
sequence of radius values (in pixels).  The \code{apMag} vector
contains the measured magnitude at any of these radius: this is the
curve-of-growth trend.  The remaining entries summaries the
relationship: \code{refRadius} is the global reference radius used for
this image; \code{maxRadius} is the outermost radius at which the
curve of growth was measured; \code{fitMag} is the fitted PSF model
magnitude integrated to infinity; \code{apRef} is the aperture
magnitude at the reference radius; \code{apLoss} is the difference
between the aperture magnitude at the reference radius and the fitted
model magnitude.  A few related functions are specified to interact
with the curve of growth:

\begin{prototype}
pmGrowthCurve *pmGrowthCurveAlloc (psF32 minRadius, psF32 maxRadius, psF32 dRadius);
\end{prototype}
This function allocates a \code{pmGrowthCurve} structure and fills in
the \code{radius} vector (see psLib SDRS \code{psVectorCreate}).  It
does {\em not} allocate the \code{apMag} vector.

\begin{prototype}
psF32 pmGrowthCurveCorrect (pmGrowthCurve *growth, psF32 radius);
\end{prototype}
This function accepts a \code{growth} curve structure and returns the
correction between the specified radius and the reference radius
($apMag(refRadius) - apMag(radius)$).

The following two functions are used to search the growth curve to the
corresponding radius entry:
\begin{prototype}
int psVectorBracket (psVector *index, psF32 key, bool above);
psF32 psVectorInterpolate (psVector *index, psVector *value, psF32 key);
\end{prototype}

\subsubsection{Aperture Trends}

With PSF model fitting, there is always some discrepancy between the
model of the PSF and the actual PSF.  As a result, the measured flux
from the model will not represent exactly the flux of the source.  It
is necessary to measure the correction between the model and the
actual source flux.  One way to perform this measurement is to compare
the model flux with the flux measured for bright stars within a fixed
aperture.  The quantity to be measured is $dA = m_{\rm aperture} -
m_{\rm fit}$.  In practice, $dA$ exhibits variations as a function of
the source position ($x,y$) and the source flux.  The variations as a
function of source position can be understood as a change in the PSF
model error as a function of position due to the changing shape of the
PSF (despite the varying PSF model, it is possible that the fitted
model yields positional variations in the residual flux).  The
variations in $dA$ as a function of magnitude can be understood as the
result of a bias in the local background measurement (for the fainter
sources) and as a result of non-linearity in the detector setting on
the bright end.  We use a 4D polynomial to represent these trends.
The first two dimensions of the polynomial represent the variation of
$dA$ as a function of $x,y$; we provide helper functions to define 1st
and 2nd order polynomials in $x,y$.  The next two dimensions are
fitted independently (no cross terms).  The first represents the
variation as a function of $r^2 / flux$, where $r$ is the aperture
radius used to measure $dA$; this is the scaling of a magnitude error
in the presence of a constant error in the sky level.  The last
dimension represents the variation of $dA$ as a function of the
stellar flux.

The following forms of the aperture correction model may be selected
by the user:
\begin{datatype}
typedef enum {
    PM_PSF_NONE,
    PM_PSF_CONSTANT,
    PM_PSF_SKYBIAS,
    PM_PSF_SKYSAT,
    PM_PSF_XY_LIN,
    PM_PSF_XY_QUAD,
    PM_PSF_SKY_XY_LIN,
    PM_PSF_SKYSAT_XY_LIN,
    PM_PSF_ALL
} pmPSF_ApTrendOptions;
\end{datatype}

The following utility function sets the aperture correction model
coefficient masks to select the specific desired coefficients:
\begin{prototype}
bool pmPSF_MaskApTrend (pmPSF *psf, pmPSF_ApTrendOptions option);
\end{prototype}

\subsubsection{pmPSF, pmPSFtry, and PSF model} 

It is useful to generate a model to define the point-spread-function
which describes the flux distribution for unresolved sources in an
image.  In general, the PSF varies with position in the image.  We
allow any of the source models defined for the \code{pmModel} to
represent the PSF.  For a given source model, the 2D spatial variation
of all of the source parameters, except the first four PSF-independent
parameters, are represented as polynomial, stored in a \code{psArray}.
The structure also contains the aperture correction model
(\code{ApTrend}) and the curve-of-growth model (\code{growth}).  The
additional elements are: \code{ApResid}, the constant term in the
aperture correction model; \code{dApResid}, the residual scatter for
bright sources ($S/N > 100$) after applying the aperture correction;
\code{skyBias}, the measured average bias in the sky measurement;
\code{skySat}, the scaling of the flux-dependent portion of the
correction.

The other elements of the structure define the quality of the PSF
determination.  

\begin{datatype}
typedef struct {
    pmModelType type;                   ///< PSF Model in use
    psArray *params;                    ///< Model parameters (psPolynomial2D)
    psPolynomial4D *ApTrend;            ///< ApResid vs (x,y,rflux) (rflux = ten(0.4*mInst)
    pmGrowthCurve *growth;              ///< apMag vs Radius
    float ApResid;                      ///< ???
    float dApResid;                     ///< ???
    float skyBias;                      ///< ???
    float skySat;                       ///< ???
    float chisq;                        ///< PSF goodness statistic
    int nPSFstars;                      ///< number of stars used to measure PSF
    int nApResid;                       ///< number of stars used to measure ApResid
} pmPSF;
\end{datatype}

\begin{prototype}
pmModel *pmModelFromPSF (pmModel *model, pmPSF *psf);
\end{prototype}
This function constructs a \code{pmModel} instance based on the
\code{pmPSF} description of the PSF.  The input is a \code{pmModel}
with at least the values of the centroid coordinates (possibly
normalization if this is needed) defined.  The values of the
PSF-dependent parameters are specified for the specific realization
based on the coordinates of the object.  

\begin{prototype}
bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask);
\end{prototype}
This function takes a collection of \code{pmModel} fitted models from
across a single image and builds a \code{pmPSF} representation of the
PSF.  The input array of model fits may consist of entries to be
ignored (noted by a non-zero \code{mask} entry).  The analysis of the
models fits a 2D polynomial for each parameter to the collection of
model parameters as a function of position (and normalization?).  In
this process, some of the input models may be marked as outliers and
excluded from the fit.  These elements will be marked with a specific
mask value (1 == \code{PSFTRY_MASK_OUTLIER}).  

We definet he following two functions to convert the PSF model
parameters into a collection of elements on a metadata structure, and
vice versa.  These can be used to read and write PSFs to a file and or
a database.
\begin{prototype}
psMetadata *pmPSFtoMD (psMetadata *metadata, pmPSF *psf);
pmPSF *pmPSFfromMD (psMetadata *metadata);
\end{prototype}

We have the capability to test several different model functions in an
attempt to build an accurate PSF for an image.  The complete set of
data needed to build and test as specific PSF model is carried by the
\code{pmPSFtry} structure:
\begin{datatype}
typedef struct {
    pmModelType modelType;
    pmPSF      *psf;
    psArray    *sources;      // pointers to the original sources
    psArray    *modelEXT;     // model fits, floating parameters 
    psArray    *modelPSF;     // model fits, PSF parameters
    psVector   *mask;
    psVector   *metric;
    psVector   *fitMag;
} pmPSFtry;
\end{datatype}
This structure contains a pointer to the collection of \code{sources}
which will be used to test the PSF model form.  It lists the
\code{pmModelType type} of model being tests, and contains an element
to store the resulting \code{psf} representation.  In addition, this
structure carries the complete collection of FLT (floating parameter)
and PSF (fixed parameter) model fits to each of the sources
\code{modelFLT} and \code{modelPSF}.  It also contains a mask which is
set by the model fitting and psf fitting steps.  For each model, the
value of the quality metric is stored in the vector \code{metric} and
the fitted instrumental magnitude is stored in \code{fitMag}.  The
quality metric for the PSF model is the aperture magnitude minus the
fitted magnitude for each source.  

This collection of aperture residuals is examined in the analysis
process, and a linear trend of the residual with the inverse object
flux (ie, $10^{0.4*mag}$) is fitted.  The result of this fit is a
measured sky bias (systematic error in the sky measured by the fits),
an effective infinite-magnitude aperture correction (\code{ApResid}),
and the scatter of the aperture correction for the ensemble of PSF
stars (\code{dApResid}).  The ultimate metric to intercompare multiple
types of PSF models is the value of the aperture correction scatter.

The following functions are used to try out a single PSF model.
\begin{prototype}
pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS);
\end{prototype}
This function takes the input collection of sources and performs a
complete analysis to determine a PSF model of the given type
(specified by model name).  The result is a \code{pmPSFtry} with the
results of the analysis.

\begin{prototype}
bool pmPSFtryMetric (pmPSFtry *try, float RADIUS);
\end{prototype}
This function is used to measure the PSF model metric for the set of
results contained in the \code{pmPSFtry} structure.

The following datatype defines the masks used by the \code{pmPSFtry}
analysis to identify sources which should or should not be included in
the analysis.
\begin{datatype}
enum {
    PSFTRY_MASK_CLEAR    = 0x00,
    PSFTRY_MASK_OUTLIER  = 0x01, // 1: outlier in psf polynomial fit (provided by psPolynomials)
    PSFTRY_MASK_EXT_FAIL = 0x02, // 2: ext model failed to converge 
    PSFTRY_MASK_PSF_FAIL = 0x04, // 3: psf model failed to converge 
    PSFTRY_MASK_BAD_PHOT = 0x08, // 4: invalid source photometry           
    PSFTRY_MASK_ALL      = 0x0f,
} pmPSFtryMaskValues;
\end{datatype}


\begin{datatype}
typedef enum {
  PM_CONTOUR_CRUDE
} pmContourType; 
\end{datatype}

Allocators for the above structures are defined as follows:
\begin{prototype}
pmSource   *pmSourceAlloc ();
pmPeak     *pmPeakAlloc (int x, int y, float counts, psPeakType class);
pmMoments  *pmMomentsAlloc ();
pmModel    *pmModelAlloc (pmModelType type);
\end{prototype}

\subsection{Basic Object Detection APIs}

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{prototype}
psVector *pmFindVectorPeaks(const psVector *vector, float threshold);
\end{prototype}

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 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 preceding 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{prototype}
psArray *pmFindImagePeaks(const psImage *image, float threshold);
\end{prototype}

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 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 an array of \code{pmPeak} entries.  The resulting set
of peaks should be considered a starting point, not an unambiguous
sample of the only real peaks.  If the input image is a subimage, the
output peak coordinates should be in the {\em parent} coordinate
frame.

\begin{prototype}
psArray *pmPeaksSubset(psArray *peaks, float maxvalue, const psRegion valid);
\end{prototype}

Create a new peaks array, removing certain types of peaks from the
input array 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.  The result of the
function is a new array with a reduced number of peaks.

\begin{prototype}
bool pmSourceDefinePixels(pmSource *mySource, 
                          pmReadout *readout,
                          psF32 x, 
                          psF32 y,
                          psF32 Radius)

bool pmSourceRedefinePixels(pmSource *mySource, 
                            pmReadout *readout,
                            psF32 x, 
                            psF32 y,
                            psF32 Radius)
\end{prototype}

The first form defines \code{psImage} subarrays (pixel, weight, and
mask) for the source located at coordinates \code{x,y} on the image
set defined by \code{readout} (in parent coords).  The pixels defined
by this operation consist of a square window (of full width $2 Radius
+ 1$) centered on the pixel which contains the given coordinate, in
the frame of the readout.  The window is defined to have limits which
are valid within the boundary of the \code{readout} image, thus if the
radius would fall outside the image pixels, the subimage is truncated
to only consist of valid pixels.  If \code{readout->mask} or
\code{readout->weight} are not \code{NULL}, matching subimages are
defined for those images as well.  This function fails if no valid
pixels can be defined (x or y less than Radius, for example).  This
function should be used to define a region of interest around a
source, including both source and sky pixels.  The second form accepts
an existing source and redefines the pixels if the requested radius
encompasses more pixels than the existing images.

\begin{prototype}
pmSource *pmSourceLocalSky(pmSource *source,
                           psStatsOptions statsOptions,
                           float Radius)
\end{prototype}

Measure the local sky in the vicinity of the given \code{source}.  The
\code{Radius} defines the square aperture in which the moments will be
measured.  This function assumes the source pixels have been defined,
and that the value of \code{Radius} here is smaller than the value of
\code{Radius} used to define the pixels.  The annular region not
contained within the radius defined here is used to measure the local
background in the vicinity of the source.  The local background
measurement uses the specified statistic passed in via the
\code{statsOptions} entry.  This function allocates the
\code{pmMoments} structure.  The resulting sky is used to set the
value of the \code{pmMoments.sky} element of the provided
\code{pmSource} structure.  

\begin{prototype}
bool pmSourceMoments(pmSource *source, float radius);
\end{prototype}

Measure source moments for the given \code{source}, using the value of
\code{source.moments.sky} provided as the local background value and
the peak coordinates as the initial source location.  The resulting
moment values are applied to the \code{source.moments} entry, and the
source is returned.  The moments are measured within the given
circular radius of the \code{source.peak} coordinates.  The return
value indicates the success (TRUE) of the operation.  This function
also measures the approximate signal-to-noise ratio of the source
(\code{source.SN}) based on the total number of source counts divided
by the square-root of the total source variance, as determined from
the weight image.

\begin{prototype}
pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata);
\end{prototype}

We use the source moments to make an initial, approximate source
classification, and as part of the information needed to build a PSF
model for the image.  As long as the PSF shape does not vary
excessively across the image, the sources which are represented by a
PSF (the start) will have very similar second moments.  The function
\code{pmSourcePSFClump} searches a collection of \code{sources} with
measured moments for a group with moments which are all very similar.
The function returns a \code{pmPSFClump} structure, representing the
centroid and size of the clump in the $\sigma_x$, $\sigma_y$
second-moment plane.  

The goal is to identify and characterize the stellar clump within the
$\sigma_x, \sigma_y$ plane.  To do this, an image is constructed to
represent this plane.  The units of $\sigma_x$ and $\sigma_y$ are in
image pixels.  A pixel in this analysis image represents 0.1 pixels in
the input image.  The dimensions of the image need only be 10 pixels.
The peak pixel in this image (above a threshold of half of the image
maximum) is found.  The coordinates of this peak pixel represent the
2D mode of the $\sigma_x, \sigma_y$ distribution.  The sources with
$\sigma_x, \sigma_y$ within 0.2 pixels of this value are then used to
calculate the median and standard deviation of the $\sigma_x,
\sigma_y$ values.  These resulting values are returned via the
\code{pmPSFClump} structure.

The return value indicates the success (TRUE) of the operation.

\tbd{limit the S/N of the candidate sources (part of Metadata)?} 

\tbd{save the clump parameters on the Metadata} 

\begin{prototype}
bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump)
\end{prototype}

Based on the specified data values, make a guess at the source
classification.  The sources are provides as a \code{psArray} of
\code{pmSource} entries.  Definable parameters needed to make the
classification are provided to the routine with the \code{psMetadata}
structure.  The rules below refer to values which can be extracted
from the metadata using the given keywords.  Except as noted, the data
type for these parameters are \code{psF32}.

The following rules are used to make the classification.  The number
of saturated pixels are counted, based on the mask having the
\code{PSPHOT_MASK_SATURATED} bit set.  Sources which are greater than
1$\sigma$ larger than the \code{pmPSFClump} center in both dimensions
and which have more than a single saturated pixel are identified as
being a likely saturated star (\code{type = PM_SOURCE_STAR, mode =
PM_SOURCE_SATSTAR}).  Sources which are not so large but which have
multiple saturated pixels are identified as saturated regions, ie
bleed trails or hot columns (\code{type = PM_SOURCE_SATURATED}).

Sources with 
\[ \sigma_x < 0.05 \]
or
\[ \sigma_y < 0.05\]
should be identified as type \code{PM_SOURCE_DEFECT} (likely cosmic
ray pixel).

Sources with
\[ \sigma_x > \mbox{CLUMP}_{x} + 3\mbox{CLUMP}_{dx}\]
and 
\[ \sigma_y > \mbox{CLUMP}_{y} + 3\mbox{CLUMP}_{dy}\]
should be identified as type \code{PM_SOURCE_EXTENDED}.  

All other sources should be identified as type \code{PM_SOURCE_STAR}.
Of these sources, the mode should be set to \code{PM_SOURCE_PSFSTAR}
for any sources with $SN$ greater than \code{PSF_SN_LIM} which are
within 1.5$\sigma$ of the PSF clump center.  These sources are used to
determine a guess at the shape of the PSF, based on the collection of
$\sigma_x$ and $\sigma_y$ values.

\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{prototype}
bool pmSourceModelGuess(pmSource *source, const psImage *image, pmModelType model);
\end{prototype}

Convert available data to an initial guess for the given model.  This
function allocates a \code{pmModel} entry for the \code{pmSource}
structure based on the provided model selection.  The method of
defining the model parameter guesses are determined by using
\code{pmModelGuessFunc_GetFunction} to determine the guess function
for the model of interest.  The returned function is called and the
guess values are used to set the model parameters.  The function
returns \code{TRUE} on success or \code{FALSE} on failure.

\begin{prototype}
psArray *pmSourceContour(const pmSource *source, const psImage *image, float level, pmContourType type);
\end{prototype}

Find points in a contour for the given source at the given level.  If
\code{type} is \code{PM_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 returned 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 contour types may be specified in the future for more refined contours}

\begin{prototype}
bool pmSourceFitModel(pmSource *source, psImage *image);
\end{prototype}

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{prototype}
bool pmModelFitStatus (pmModel *model);
\end{prototype}

This function wraps the call to the model-specific function returned
by \code{pmModelFitStatusFunc_GetFunction}.  The model-specific
function examines the model parameters, parameter errors, Chisq, S/N,
and other parameters available from \code{model} to decide if the
particular fit was successful or not.

\begin{prototype}
bool pmSourceAddModel(psImage *image, pmSource *source, bool center, bool sky);
bool pmSourceSubModel(psImage *image, pmSource *source, bool center, bool sky);
\end{prototype}

Add or subtract the given source model flux to/from the provided
image.  The boolean option \code{center} selects if the source is
re-centered to the image center or if it is placed at its centroid
location.  The boolean option \code{sky} selects if the background sky
is applied (\code{TRUE}) or not.  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{prototype}
bool pmSourcePhotometry (float *fitMag,  // integrated fit magnitude
                         float *obsMag,  // aperture flux magnitude
                         pmModel *model, // model used for photometry
                         psImage *image, // image pixels to be used
                         psImage *mask   // mask of pixels to ignore
);
\end{prototype}

The function returns both the magnitude of the fit, defined as $-2.5
\log{\rm flux}$, where the flux is integrated under the model,
theoretically from a radius of 0 to infinity.  In practice, we
integrate the model beyond $50 \sigma$.  The aperture magnitude is
defined as $-2.5 \log{\rm flux}$, where the flux is summed for all
pixels which are not excluded by the aperture mask.  The model flux is
calculated by calling the model-specific function provided by
\code{pmModelFlux_GetFunction}.

\begin{prototype}
int pmSourceDophotType (pmSource *source);
\end{prototype}
This function converts the source classification into the closest
available approximation to the Dophot classification scheme.  The
following list gives the correspondence:
\begin{verbatim}
PM_SOURCE_DEFECT:       8
PM_SOURCE_SATURATED:    8
PM_SOURCE_SATSTAR:      10
PM_SOURCE_PSFSTAR:      1
PM_SOURCE_GOODSTAR:     1
PM_SOURCE_POOR_FIT_PSF: 7
PM_SOURCE_FAIL_FIT_PSF: 4
PM_SOURCE_FAINTSTAR:    4
PM_SOURCE_GALAXY:       2
PM_SOURCE_FAINT_GALAXY: 2
PM_SOURCE_DROP_GALAXY:  2
PM_SOURCE_FAIL_FIT_GAL: 2
PM_SOURCE_POOR_FIT_GAL: 2
PM_SOURCE_OTHER:        ?
\end{verbatim}

\begin{prototype}
int pmSourceSextractType (pmSource *source);
\end{prototype}
This function converts the source classification into the closest
available approximation to the Sextractor classification scheme.
\tbd{the correspondence is not yet defined}.

\subsection{Object List Input/Output}

We support several object catalog formats.  Some of these mimic the
formats used by the Elixir system to support testing with existing
data and software.  Some of these are for use by the Pan-STARRS
project for testing.

\subsubsection{OBJ Format}

This format is produced by versions of DoPhot and is used by the
Elixir system as an intermediate output data product.  The objects are
written to a text file with fixed line-length and with fixed column
positions.  The file has no header associated with it.  This is only
an output format, and should be used just for testing and comparison
with the Elixir tools.

\subsubsection{SX Format}

This format is produced by versions of Sextractor and is used by the
Elixir system as an intermediate output data product.  The objects are
written to a text file with fixed line-length and with fixed column
positions.  The file has no header associated with it.  This is only
an output format, and should be used just for testing and comparison
with the Elixir tools.  The SX and OBJ formats are similar, but use a
somewhat different definition of the columns.

\subsubsection{CMP Format}

This format is used extensively by the Elixir system, and many data
files are available in this format.  The format is a pseudo-FITS
format, consisting of a FITS header (with NAXIS=2) and a text data
segment with fixed line length.  The CMP files are always in SPLIT
format in the sense that each object table is a single file. 

\subsubsection{CMF Format}

This format is a true FITS table format.  The object data is stored
for each readout in a separate extension.  In addition, the Cell
headers are stored in their own extensions (with NAXIS=0).  In SPLIT
format, the Cell header is the PHU header.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Image Combination}

\tbd{PAP to code and update.}

The image combination for \PS{} will employ an iterative approach, in
order to identify cosmic rays.  The first pass involves transforming
and combining the input images, and noting pixels which are apparently
deviant.  These pixels are examined in further detail, before a subset
of them are declared to be bad, whereupon these pixels are
re-transformed, and the images are combined properly.  Here we
introduce two functions which will perform the combination and
examination steps.  Prototype code exists for each of these functions.
\tbd{For further details, see the document about image combination for
\PS{}.}

\subsection{Combining images}

\begin{prototype}
psImage *pmCombineImages(psImage *combined, // Combined image
                         psArray **questionablePixels, // Array of rejection masks
                         const psArray *images, // Array of input images
                         const psArray *errors, // Array of input error images
                         const psArray *masks,// Array of input masks
                         unsigned int maskVal, // Mask value
                         const psPixels *pixels, // Pixels to combine
                         int numIter,   // Number of rejection iterations
                         float sigmaClip, // Number of standard deviations at which to reject
                         const psStats *stats // Statistics to use in the combination
                         );
\end{prototype}

\code{pmCombineImages} shall combine the input \code{images},
returning the \code{combined} image and a list of
\code{questionablePixels} in each input image.  The array of error
images, \code{errors}, shall be used to calculate the value in the
combined image and the list of questionable pixels, if
non-\code{NULL}.  Pixels whose corresponding value in the array of
mask images, \code{masks}, matches \code{maskVal} shall be masked from
the combination.  The \code{images}, \code{errors} and \code{masks}
arrays, if non-\code{NULL}, shall all carry the same number of images;
otherwise the function shall generate an error and return \code{NULL}.
The sizes of all images in the \code{images}, \code{errors} and
\code{masks} arrays shall be identical; otherwise the function shall
generate an error and return \code{NULL}.

If \code{pixels} is non-\code{NULL}, only those pixels specified shall
be combined.  The combination consists of \code{numIter} iterations in
which a stack of pixels is combined using the specified \code{stats}.
In each iteration, questionable pixels are identified as lying more
than \code{sigmaClip} standard deviations from the combined value;
these pixels are excluded from the stack for the next iteration.  The
value for the combined image is that produced by the \textit{first}
iteration (i.e., with no pixels excluded except those which have their
corresponding mask match the \code{maskVal}); this allows subsequent
calls to the function to only act on a small fraction of the pixels,
since questionable pixels identified in the first call of the function
will be properly rejected at a later point (see the example, below).

In the event that \code{images} or \code{stats} are \code{NULL}, the
function shall generate an error and return \code{NULL}.

\subsection{Rejecting pixels}

\begin{prototype}
psArray *pmRejectPixels(const psArray *images, // Array of input images
                        const psArray *masks, // Array of masks for input images
                        const psArray *pixels, // These are the pixels which were rejected in the combination
                        const psArray *inToOut, // Transformations from input to output system
                        const psArray *outToIn, // Transformations from output to input system
                        float rejThreshold, // Rejection threshold
                        float gradLimit // Gradient limit
                        );
\end{prototype}

\tbd{This algorithm will change: an addition will be made to avoid
masking pixels in the wings of a star when combining images taken in
different seeing, and the gradient limit criteria will be changed.}

\code{pmRejectPixels} inspects those questionable \code{pixels}
identified by \code{pmCombineImages} to determine if they are truly
discrepant.  This inspection is performed in the coordinate frame of
the detector, where the pixels haven't been smeared by transformation.
Two tests are applied to each of the \code{images}:
\begin{enumerate}
\item The list of questionable pixels for an image is converted to an
  image which is transformed back to the coordinate frame of the
  detector.  Those pixels in the detector frame which have a value
  exceeding \code{rejThreshold} are suspected cosmic rays and
  subjected to the next test.  Depending on the value of the
  \code{rejThreshold}, this test basically amounts to demanding that
  questionable pixels neighbor each other in the transformed image.
\item The cores of point sources may mimic a cosmic ray, especially in
  under-sampled images.  To minimize flagging stars as cosmic rays, we
  determine the gradient around the pixel of interest; if the gradient
  is large, then the pixel is likely the core of a point source.  In
  order to reliably measure the gradient in the presence of a
  suspected cosmic ray, we use the companion images --- the gradient
  is the mean gradient at the corresponding position on the other
  images.  In order to calculate the corresponding positions, the
  \code{inToOut} and \code{outToIn} transformations are required.  If
  the gradient is less than \code{gradLimit}, then the pixel is
  identified as a cosmic ray.
\end{enumerate}

The function shall return an array of \code{psPixels}, one for each of
the input \code{images}, containing pixels that have been identified
as cosmic rays according to the above criteria.

If any of the input pointers are \code{NULL}, then the function shall
generate an error and return \code{NULL}.

\subsection{Example}

Here is an example of what the image combination routine looks like,
demonstrating how the various pieces fit together.  The inputs are:
\begin{itemize}
\item \code{psArray *inputs}: Input detector images, each a
  \code{psImage} of type \code{psF32}
\item \code{psArray *inputMask}: Input mask images, each a
  \code{psImage} of type \code{psU8}
\item \code{psArray *inputsErr}: Input error images, each a
  \code{psImage} of type \code{psF32}
\item \code{psPlaneTransform *skyToDetector}: Maps from sky
  coordinates to detector coordinates, each a \code{psPlaneTransform}
\item \code{psRegion *combineRegion}: Sky coordinate pixels to combine
\item \code{int numIter}: Number of iterations in combination
\item \code{float rejThreshold}: Threshold for rejection
\item \code{float gradLimit}: Limit for gradient
\end{itemize}

The output is the combined image.

\begin{verbatim}
    psArray *transformed = psArrayAlloc(nImages); // Array of transformed images
    psArray *transformedErr = psArrayAlloc(nImages); // Array of transformed error images
    psArray *transformedMask = psArrayAlloc(nImages); // Array of masks for transformed images

    for (int i = 0; i < nImages; i++) {
        psPixels *blanks = NULL;        // List of blank pixels
        transformed->data[i] = psImageTransform(NULL, &blanks, inputs->data[i],
                                                inputMask->data[i], inputMaskVal, NAN, skyToDetector,
                                                combineRegion, NULL, PS_INTERPOLATE_BILINEAR);
        transformedErr->data[i] = psImageTransform(NULL, NULL, inputsErr->data[i], inputMask->data[i],
                                                   inputMaskVal, NAN, skyToDetector, combineRegion, NULL,
                                                   PS_INTERPOLATE_BILINEAR_VARIANCE);
        psImage *skyImage = transformed->data[i]; // Dereference the transformed image
        psRegion *blankRegion = psRegionAlloc(0, 0, skyImage->numCols, skyImage->numRows); // Size of
                                                                                           // transformed
                                                                                           // image
        transformedMask->data[i] = psPixelsToMask(NULL, blanks, *blankRegion, PS_MASK_BLANK);
        psFree(blankRegion);
        psFree(blanks);
    }

    psArray *rejected = NULL;           // Array of rejected pixel lists
    psStats *combineStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistic to use in doing the combination
    psImage *combined = pmCombineImages(NULL, &rejected, transformed, transformedErr, transformedMask, 0,
                                        NULL, numIter, sigmaClip, combineStats); // Combined image
    psArray *bad = pmRejectPixels(inputs, rejected, NULL, skyToDetector, rejThreshold, gradLimit); // Bad pix
    psPixels *combinePixels = NULL;     // Pixels to combine
    for (int i = 0; i < nImages; i++) {
        psPixels *badSource = psPixelsTransform(NULL, bad->data[i], skyToDetector); // Bad pixels on the input
        psImage *badMask = psPixelsToMask(NULL, badSource, PS_MASK_COSMICRAY); // Mask image for the input
        (void)psBinaryOp(inputMask->data[i], inputMask->data[i], "|", badMask); // Put CRs into original mask
        psFree(badSource);
        psFree(badMask);

        combinePixels = psPixelsConcatenate(redo, bad->data[i]);

        // Update transformed image
        psPixels *blanks = NULL;        // List of blank pixels
        transformed->data[i] = psImageTransform(transformed->data[i], &blanks, inputs->data[i],
                                                inputMask->data[i], inputMaskVal | PS_MASK_COSMICRAY, NAN,
                                                skyToDetector, combineRegion, bad->data[i],
                                                PS_INTERPOLATE_BILINEAR);
        transformedErr->data[i] = psImageTransform(transformedErr->data[i], NULL, inputsErr->data[i],
                                                   inputMask->data[i], inputMaskVal | PS_MASK_COSMICRAY,
                                                   NAN, skyToDetector, combineRegion, bad->data[i],
                                                   PS_INTERPOLATE_BILINEAR_VARIANCE);
        psImage *skyImage = transformed->data[i]; // Dereference the transformed image
        psRegion *blankRegion = psRegionAlloc(0, 0, skyImage->numCols, skyImage->numRows); // Size of
                                                                                           // transformed
                                                                                           // image
        transformedMask->data[i] = psPixelsToMask(transformedMask->data[i], blanks, *blankRegion,
                                                  PS_MASK_BLANK);
        psFree(blankRegion);
        psFree(blanks);
    }
    psFree(bad);

    // Combine with no rejection
    combined = pmCombineImages(combined, NULL, transformed, transformedErr, transformedMask,
                               PS_MASK_BLANK, combinePixels, 0, 0.0, combineStats);
    psFree(combineStats);
    psFree(combinePixels);
    psFree(transformed);
    psFree(transformedErr);
    psFree(transformedMask);
\end{verbatim}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Image Subtraction}

\tbd{PAP to code and update.}

Image subtraction is arguably the best method of identifying faint
variable sources in images with different point-spread functions.  It
relies on fitting for a convolution kernel that minimizes the
residuals in subtracting small regions of the image.  The use of a
convolution kernel consisting of a linear combination of basis
functions allows the problem to be solved with only modest computing
power.

\subsection{The kernels}

We will allow for the use of two convolution kernels.  The first is
that employed by the popular image subtraction program,
\href{http://www2.iap.fr/users/alard/package.html}{ISIS}, consisting
of Gaussians modified by polynomials:
\begin{equation}
B_{ijk}(u,v) = e^{-(u^2 + v^2)/2\sigma_i^2} u^j v^k
\end{equation}
The second simply consists of delta functions, which we refer to as
POIS (Pan-STARRS Optimal Image Subtraction):
\begin{equation}
B_{ij}(u,v) = \delta(u - i)\ \delta(v - j)
\end{equation}
\tbd{For further details, see the document about image subtraction for
\PS{}.}  The former is widely used, while the second appears to be
equally useful and faster, though not as tried and proven.

\begin{datatype}
typedef enum {
    PM_SUBTRACTION_KERNEL_POIS,         // POIS kernel --- delta functions
    PM_SUBTRACTION_KERNEL_ISIS          // ISIS kernel --- gaussians modified by polynomials
} pmSubtractionKernelsType;
\end{datatype}

In order to simplify the book-keeping for the kernels, we will define
a \code{pmSubtractionKernels}, which keeps track of the details of the
each of the kernel basis functions:

\begin{datatype}
typedef struct {
    pmSubtractionKernelType type;       // Type of kernels --- allowing the use of multiple kernels
    int size;                           // Size of kernel in x and y
    int spatialOrder;                   // Maximum order of spatial variations
    psVector *u, *v;                    // Offset (for POIS) or polynomial order (for ISIS)
    psVector *sigma;                    // Width of Gaussian (for ISIS)
    psVector *xOrder, *yOrder;          // Spatial polynomial order (for all)
    int subIndex;                       // Index of kernel to be subtracted to maintain flux conservation
    psArray *preCalc;                   // Array of images containing pre-calculated kernel (to
                                        // accelerate ISIS; don't use for POIS)
} pmSubtractionKernels;
\end{datatype}

This structure caters for both choices of kernel type.  For a POIS
kernel, the \code{u} and \code{v} vectors shall be set to the
coordinates for the delta functions for the corresponding kernel.  For
an ISIS kernel, the \code{sigma} vector shall be set to the Gaussian
widths and the \code{u} and \code{v} vectors shall be set to the
orders of the modifying polynomials for the corresponding kernel.  For
both choices of kernel, the \code{xOrder} and \code{yOrder} vectors
specify the order of the spatial variation.

In order to maintain flux conservation when the kernel is spatially
variable, we need to treat one kernel in the set differently.  The
convolutions for this kernel, identified by the \code{subIndex}, are
calculated in the usual way, while all others have the \code{subIndex}
kernel subtracted from them.  For details, see the
\href{http://www.edpsciences.org/journal/index.cfm?v_url=aas/full/2000/11/ds8706/ds8706.html}{paper
by Alard (2000, A\&AS, 144, 363)}.

Since the ISIS kernels are continuous functions, it is worth
pre-calculating them instead of calculating them each time they are
required.  The \code{preCalc} array, consisting of \code{psImage}s is
provided for this purpose.

The \code{pmSubtractionKernels} are generated by the following functions:

\begin{prototype}
pmSubtractionKernels *pmSubtractionKernelsAllocPOIS(int size, int spatialOrder);
pmSubtractionKernels *pmSubtractionKernelsAllocISIS(const psVector *sigmas, const psVector *orders,
                                                    int size, int spatialOrder);
\end{prototype}

\code{pmSubtractionKernelsAllocPOIS} shall generate the
\code{pmSubtractionKernels} suitable for the POIS kernel basis set.
This involves setting the \code{u}, \code{v}, \code{xOrder} and
\code{yOrder} to the appropriate values.  \code{size} is the half-size
of the kernel, and \code{spatialOrder} is the maximum spatial order
(the spatial variation is $x^i y^j$ with $i+j <$ \code{spatialOrder}).
The \code{subIndex} is set to the kernel which has \code{u = 0},
\code{v = 0}, \code{xOrder = 0} and \code{yOrder = 0}.  There should
be \code{(2 * size + 1) * (2 * size + 1) * (spatialOrder + 1) *
(spatialOrder + 2) / 2} kernels.

\code{pmSubtractionKernelsAllocISIS} shall generate the
\code{pmSubtractionKernels} suitable for the ISIS kernel basis set.
This involves setting the \code{sigma}, \code{u}, \code{v},
\code{xOrder} and \code{yOrder} to the appropriate values, as well as
generating the \code{preCalc} images.  Note that the \code{sigma}
vector contained within the \code{pmSubtractionKernels} is not the
same as the input \code{sigmas} vector, but contains repeated entries.
\code{size} is the half-size of the kernel, which specifies the size
of the \code{preCalc} images.  The \code{spatialOrder} is the maximum
spatial order (the spatial variation is $x^i y^j$ with $i+j <$
\code{spatialOrder}).  The \code{subIndex} is set to the kernel which
has \code{u = 0}, \code{v = 0}, \code{xOrder = 0} and \code{yOrder =
0}, for the first of the Gaussian widths in the \code{sigmas} vector.

\subsection{Stamps}

Sub-regions on an image which are used to derive the best-fit
convolution kernel are referred to as ``stamps''.

\begin{datatype}
typedef struct {
    int x, y;                           // Position
    psImage *matrix;                    // Associated matrix
    psVector *vector;                   // Associated vector
    pmStampStatus status;               // Status of stamp
} pmStamp;
\end{datatype}

A stamp is the region around a central pixel, \code{x,y}.  The
\code{matrix} and \code{vector} are generated in the process of
solving for the best-fit convolution kernel; each of these will likely
be of type \code{psF64} in order to maintain the best possible
precision (we will be summing squares).  In order to allow us to throw
out stamps without having to laboriously recompute the total
least-squares matrix and vector, we use a separate matrix and vector
for each stamp.

To allow iteration on the choice of stamps, a stamp contains a
\code{status}, an enumerated type:

\begin{datatype}
typedef enum {
    PM_STAMP_USED,                      // Use this stamp
    PM_STAMP_REJECTED,                  // This stamp has been rejected
    PM_STAMP_RECALC,                    // Having been reset, this stamp needs to be recalculated
    PM_STAMP_NONE                       // No stamp in this region
} pmStampStatus;
\end{datatype}

\begin{prototype}
psArray *pmSubtractionFindStamps(psArray *stamps, // Output stamps, or NULL
                                 const psImage *image, // Image for which to find stamps
                                 const psImage *mask, // Mask
                                 unsigned int maskVal, // Value for mask
                                 float threshold, // Threshold for stamps in the image
                                 int xNum, int yNum, // Number of stamps in x and y
                                 int border // Border around image to ignore (should be size of kernel)
                                 );
\end{prototype}

\code{pmSubtractionFindStamps} returns an array of stamps on the
\code{image} suitable for use in calculating the best-fit convolution
kernel.  Except for a \code{border} all the way around, the
\code{image} is broken into \code{xNum} $\times$ \code{yNum}
rectangles; there will be a stamp within each rectangle.  If
\code{stamps} is non-\code{NULL}, then the function shall only attempt
to identify a new stamp in a particular rectangle if the corresponding
stamp \code{status} is \code{PM_STAMP_REJECTED}.

A stamp shall be recognized as the pixel with the greatest value that
does not have the corresponding pixel in the \code{mask} matching
\code{maskVal}.  If the value of the this pixel does not exceed
\code{threshold}, then the stamp \code{status} shall be marked as
\code{PM_STAMP_NONE}, which means that the stamp will be ignored in
future iterations.  If a legitimate stamp is found within the region,
then its status shall be changed to \code{PM_STAMP_RECALC}.


\subsection{Solving for the kernel}

Calculating the best-fit convolution kernel requires solving a matrix
equation, the elements of which are obtained by applying the kernel
basis functions to the stamps.  The final matrix and vector are the
sum of the matrices and vectors obtained for each of the individual
stamps.

\begin{prototype}
bool pmSubtractionCalculateEquation(psArray *stamps, // The stamps for which to calculate the equation
                                    const psImage *reference, // Reference image
                                    const psImage *input, // Input image
                                    const psSubtractionKernels *kernels, // The kernel basis functions
                                    int footprint // Half-size of region over which to calculate equation
                                    );
\end{prototype}

\code{pmSubtractionCalculateEquation} shall calculate the
\code{matrix} and \code{vector} for each of the \code{stamps} which
have \code{status} set to \code{PM_STAMP_RECALC}.  The calculation is
made over a region with a half size of \code{footprint} on the
\code{reference} and \code{input} images, using each of the
\code{kernels}.  In the event that any of the input pointers are
\code{NULL}, the function shall generate an error and return
\code{false}; otherwise, the function shall return \code{true}.

The vector is:
\begin{equation}
v_i = \sum_{x,y} I(x,y) [ R(x,y) \otimes B_i(u,v) ] / \sigma(x,y)^2
\end{equation}
and the matrix is:
\begin{equation}
M_{ij} = \sum_{x,y} \left[ R(x,y) \otimes B_i(u,v) \right] \  \left[ R(x,y) \otimes B_j(u,v) \right] / \sigma(x,y)^2
\end{equation}
where $I(x,y)$ is the input image, $R(x,y)$ is the reference image,
$B_i(u,v)$ is the $i$-th kernel basis function, $\otimes$ denotes
convolution, $\sigma(x,y) = R(x,y)^{1/2}$ is an estimate of the error,
and the sum over $x,y$ indicates summing over the stamp regions.

In addition to the each of the \code{kernels}, an additional parameter
for which we must solve is the difference in the background level
between the \code{reference} and \code{input} images.  The appropriate
term shall be added to the \code{matrix} and \code{vector}.

In order to maintain flux conservation when the kernel is spatially
variable, for each of the kernel basis functions apart from the first,
the kernel actually employed shall be the first kernel function
subtracted from the original kernel function.

Having calculated the matrix equation for a stamp, its \code{status}
is set to \code{PM_STAMP_USED}.

Since this step is one of the major rate-limiting factors in image
subtraction, care should be taken with optimization.

\begin{prototype}
psVector *pmSubtractionSolveEquation(psVector *solution,        // Solution vector, or NULL
                                     const psArray *stamps // Array of stamps
                                     );
\end{prototype}

\code{pmSubtractionSolveEquation} shall solve the matrix equation
provided by each of the \code{stamps}, returning the \code{solution}
vector.  This involves summing the \code{matrix} and \code{vector} of
each of the stamps which have \code{status} set to
\code{PM_STAMP_USED}, and multiplying the inverse of the matrix by the
\code{vector}.  If the \code{solution} is \code{NULL}, then the
function shall allocate and return a new vector; otherwise, the
\code{solution} vector shall be modified in-place.  If \code{stamps}
is \code{NULL}, then the function shall generate an error and return
\code{NULL}.  The type of the \code{solution} vector should be
\code{psF64}, since the matrix equation involves summing squares.


\subsection{Rejection of stamps}

\begin{prototype}
bool pmSubtractionRejectStamps(psArray *stamps, // Array of stamps to check for rejection
                               psImage *mask, // Mask image
                               unsigned int badStampMaskVal, // Value to use in mask for bad stamp
                               int footprint, // Region to mask if stamp is bad
                               float sigmaRej, // Number of RMS deviations above zero at which to reject
                               const psImage *refImage, // Reference image
                               const psImage *inImage, // Input image
                               const psVector *solution, // Solution vector
                               const pmSubtractionKernels *kernels // Array of kernel parameters
                               );
\end{prototype}

\code{pmSubtractionRejectStamps} shall apply the \code{solution} to
the \code{stamps}, rejecting stamps for which the mean square
residuals exceed \code{sigmaRej} RMS deviations from zero.
\code{stamps} which are rejected have their \code{status} set to
\code{PM_STAMP_REJECTED}, and have pixels within \code{footprint} of
the corresponding position in the \code{mask} set to
\code{badStampMaskVal} so they will not be used again.

The deviations are calculated through extracting the stamps from the
\code{refImage} and \code{inImage}, convolving the reference stamp by
the best-fit kernel (derived from the \code{solutions} vector and the
\code{kernels}), subtracting and then dividing by the stamp from the
input image, and then squaring to obtain the mean square residual.

\subsection{Visualization of kernel}

Having solved for the best-fit kernel, it is often useful to visualize
it.

\begin{prototype}
psImage *pmSubtractionKernelImage(psImage *out, const psVector *solution,
                                  const pmSubtractionKernels *kernels, float x, float y);
\end{prototype}

\code{pmSubtractionKernelImage} shall create an image of the kernel
from the \code{solution} vector and the \code{kernels}.  The relative
position (between -1 and +1) on the image at which to evaluate the
kernel (important if the kernel is spatially variable) is specified by
\code{x} and \code{y}.  If \code{out} is \code{NULL}, then the
function shall allocate a new image of sufficient size (matching the
\code{precalc} images), and return the result; otherwise, \code{out}
shall be modified in-place.


\subsection{Example}

Here is an example of what the image subtraction routine looks like,
demonstrating how the various pieces fit together.  The inputs are:
\begin{itemize}
\item \code{psImage *reference}: Reference image
\item \code{psImage *refMask}: Mask for reference image
\item \code{psImage *input}: Input image
\item \code{psImage *inMask}: Mask for input image
\item \code{unsigned int maskVal}: Value to be masked
\item \code{pmSubtractionKernelType kernelType}: Type of kernel to use
\item \code{int kernelHalfSize}: Half the kernel size (full size is \code{2*kernelHalfSize + 1})
\item \code{psVector *sigmas}: Widths for the ISIS Gaussians
\item \code{psVector *polyOrders}: Polynomial orders for ISIS Gaussians
\item \code{int spatialOrder}: Maximum spatial order for spatially variable kernel
\item \code{float stampThreshold}: Threshold for finding stamps
\item \code{int nStampsX, nStampsY}: Number of stamps in x and y
\item \code{int stampSize}: Half size of stamp footprint
\item \code{int numIter}: Number of iterations on the stamps
\item \code{float sigmaRej}: Rejection threshold for stamps
\end{itemize}

The output is the subtracted image and the corresponding mask.

\begin{verbatim}
    // Mask around bad pixels in the reference image.  There are two cases to worry about:
    // 1. Bad pixels within the kernel, which will affect the subtracted image
    // 2. Bad pixels within the stamp, which affects the calculation of the kernel
    psImage *subMask = psImageGrowMask(NULL, refMask, maskVal, kernelHalfSize, PS_MASK_NEAR_BAD);
    (void)psImageGrowMask(subMask, refMask, maskVal, stampSize, PS_MASK_BAD_STAMP);
    // Add in the mask for the input image.  Don't need to grow this, since it isn't convolved.
    (void)psBinaryOp(subMask, subMask, "|", inMask);

    // Generate kernel basis functions
    psArray *kernels = NULL;            // Array of kernel basis functions
    switch (kernelType) {
      case PM_SUBTRACTION_KERNEL_POIS:
        // Create the kernel basis functions
        kernels = pmSubtractionKernelsGeneratePOIS(kernelHalfSize, spatialOrder);
        break;
      case PM_SUBTRACTION_KERNEL_ISIS:
        kernels = pmSubtractionKernelsGenerateISIS(sigmas, polyOrders, kernelHalfSize, spatialOrder);
        break;
      default:
        barf();
    }

    psArray *stamps = NULL;             // Array of stamps
    psVector *kernelCoeffs = NULL;      // Coefficients for the kernels
    bool rejected = true;               // Did we reject a stamp in the last iteration?

    // Iterate for a solution
    for (int iter = 0; iter < numIter && rejected; iter++) {

        // Find stamps
        stamps = pmSubtractionFindStamps(stamps, reference, subMask, maskVal | PS_MASK_BAD_STAMP,
                                         stampThreshold, nStampsX, nStampsY, stampSize, kernelHalfSize);

        // Generate and solve matrix equations
        (void)pmSubtractionCalculateEquation(stamps, reference, input, kernels, stampSize);
        kernelCoeffs = pmSubtractionSolveEquation(kernelCoeffs, stamps);

        // Reject bad stamps
        rejected = pmSubtractionRejectStamps(stamps, subMask, PS_MASK_BAD_STAMP, stampSize, sigmaRej,
                                             reference, input, kernelCoeffs, kernels);
    }

    // Convolve the reference image
    psImage *referenceConvolved = pmSubtractionConvolveImage(NULL, reference, subMask, kernelCoeffs, kernels);
    // Subtract
    psImage *subtracted = (psImage*)psBinaryOp(NULL, input, "-", referenceConvolved);

    // What does the kernel look like?
    psImage *kernelImage = pmSubtractionKernelImage(NULL, kernelCoeffs, kernels, 0.0, 0.0);
    // Check/save kernel image, print statistics....

    psFree(referenceConvolved);
    psFree(stamps);
    psFree(kernels);
    psFree(kernelCoeffs);
\end{verbatim}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\appendix

\section{Basic Object Models}
\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{prototype}
float pmMinLM_Gauss2D(psVector *deriv, psVector *params, psVector *x);
\end{prototype}

This function is a two-dimensional Gaussian with an elliptical
cross-section and a constant local background.  

The initial guess for the Gaussian parameters may be taken from the
moments, peak value, and local sky.

\subsubsection{Pseudo-Gaussian}

\begin{prototype}
float pmMinLM_PseudoGauss2D(psVector *deriv, psVector *params, psVector *x);
\end{prototype}

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 initial guess for the Gaussian parameters may be taken from the
moments, peak value, and local sky.

\subsubsection{Waussian}

\begin{prototype}
float pmMinLM_Wauss2D(psVector *deriv, psVector *params, psVector *x);
\end{prototype}

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{prototype}
float pmMinLM_TwistGauss2D(psVector *deriv, psVector *params, psVector *x);
\end{prototype}

This function describes an object with power-law wings and a flattened
core, where the core has a different contour from the wings.  

The initial 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{prototype}
float pmMinLM_Sersic(psVector *deriv, psVector *params, psVector *x);
\end{prototype}

\subsubsection{Sersic with Core Galaxy Model}

\begin{prototype}
float pmMinLM_SersicCore(psVector *deriv, psVector *params, psVector *x);
\end{prototype}

\subsubsection{Pseudo Sersic Galaxy Model}

\begin{prototype}
float pmMinLM_PseudoSersic(psVector *deriv, psVector *params, psVector *x);
\end{prototype}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Revision Change Log}
\input{ChangeLogSDRS.tex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{plain}
\bibliography{panstarrs}

\end{document}

%%% All this is here for storage only --- it's not part of the document.


\subsection{Cosmic rays}

Given an input image, a choice of a particular algorithm with
corresponding parameters, \code{psPhase2MaskCRs} shall mask cosmic
rays on the input image on the basis of their morphology.  The API shall be the following:
\begin{prototype}
/** Masks Cosmic Rays on the input image on the basis of morphology. */
psReadout *psPhase2MaskCRs(psReadout *in, ///< Input image to be masked, and output
                           int algorithm, ///< Algorithm number to use
                           const void *params ///< Parameters for algorithm
                           );
\end{prototype}

This is one case in which the best choice of algorithm is not known,
and may change in the future.  For this reason, we specify the inputs
as an integer to specify the choice of algorithm in addition to a
pointer to some parameters which will be decoded on the basis of the
choice of algorithm.

Note that the input image is modified in-place.

\subsection{psPhase2MaskOpticalDefects}

Given an input image, a list of nearby stars, and a growing radius,
\code{psPhase2MaskOpticalDefects} shall mask optical defects on the
image.  The API shall be the following:
\begin{prototype}
/** Masks optical defects in the input image */
psChip *psPhase2MaskOpticalDefects(psChip *in, ///< Image to be masked (with astrometry), and output
                                   const psDlist *catalog, ///< Catalog stars nearby: a list of psObjects
                                   int grow ///< Number of pixels to grow
                                   );
\end{prototype}

\tbd{It's not clear to me how this is accomplished apart from an
optical model of the camera.  Put this one on the backburner?}


\section{Objects}

To identify and measure objects, we must measure the PSF, and then
convolve the image by this PSF and apply a threshold.

\subsection{Measure the PSF}

Given an input image, a choice of algorithm with corresponding
parameters, \code{psPhase2MeasurePSF} shall return the PSF for the
image.  The API shall be the following:
\begin{prototype}
/** Measures the PSF on the input image.  Returns the PSF */
psImage *psPhase2MeasurePSF(const psReadout *in, ///< Input image for which to measure the PSF
                            int algorithm, ///< Algorithm number to use
                            const void *params ///< Parameters for algorithm
                            );
\end{prototype}

This is another case where the algorithm is not currently clear, and
may change in the future.  For this reason, we specify the inputs to
be a choice of algorithm, and a pointer to some parameters, which are
interpreted on the basis of the algorithm choice.

\subsection{Find and measure objects}

Given an input image, the PSF of that image, and a list of flux levels
at which to threshold, \code{psPhase2GetObjects} shall return a
readout with the \code{objects} member set to a list of objects.  The
API shall be the following:
\begin{prototype}
/** Find and measure objects on the input image.  Fills in the "objects" member of the psReadout. */
/** THIS NEEDS WORK. */
psReadout *psPhase2FindObjects(psReadout *in, ///< Input image on which to find objects, and output
                               const psImage *psf, ///< PSF to use to find objects
                               const psVector *levels ///< Threshold levels (std dev.s)
                               );
\end{prototype}

Note that the input image shall be modified in-place, only insofar as
the \code{objects} member shall be set to the list of objects on that
image.

\tbd{This needs a lot more work.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Astrometry}

Given a chip, the elements of which have objects found and measured on
them, a list of catalog stars which lie on (or near) the chip, and
clipping parameters, \code{psPhase2Astrometry} shall fit an
astrometric solution.  The API shall be the following:
\begin{prototype}
/** Corrects astrometry on the input chip */
psChip *psPhase2Astrometry(psChip *in,  ///< Input chip for which to do astrometry, and output
                           const psDlist *catalog, ///< Catalog stars on the chip: a list of psObjects
                           int nClips,  ///< Number of clipping iterations
                           float clipLevel ///< Level at which to clip
                           );
\end{prototype}

Note that the input chip shall be modified in-place, only insofar as
the appropriate astrometry members shall be updated to correspond to
the fit solution (specifically, \code{in->cells[i]->cellToSky}, for
each \code{i}).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Postage stamps}

Postage stamps consist of subimages of specific objects of interest,
which are saved for more careful analysis of the pixels.

\subsection{Specification}

The postage stamps are specified by a center and size, in celestial
coordinates.  We define a \code{psPostageStampSpec} to specify the
parameters for a postage stamp:
\begin{datatype}
/** Specification of a postage stamp: location on the sky */
typedef struct {
    psSphereCoord *center;              ///< Centre of postage stamp
    psSphereCoord *size;                ///< Size of postage stamp
} psPostageStampSpec;
\end{datatype}

The \code{center} shall be specified in ICRS coordinates, which is the
\PS{} standard system.  The \code{size} shall be specified in
arcseconds on the sky.

\subsection{Extracting postage stamps}

Given an input chip, and a linked-list of regions on or near the input
chip, \code{psPhase2PostageStamps} shall output an array of subimages,
containing each of the regions.  The API shall be the following:
\begin{prototype}
/** Return postage stamps of a set of regions */
psImageArray *psPhase2PostageStamps(const psChip *in, ///< Chip from which to form postage stamps
                                    const psDlist *regions ///< Regions to postage-stampise: a list of
                                                           ///< psPostageStampSpec
                                    );
\end{prototype}

\code{regions} shall be a linked list of \code{psPostageStampSpec}s,
not all of which may correspond to legal positions on the input chip,
\code{in}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Fixed Pattern}

The fixed pattern is a correction to the general astrometric solution
formed by summing the residuals from many observations.  The intent is
to correct for higher-order distortions in the camera system on a
coarse grid (larger than individual pixels, but smaller than a single
cell).  Hence, in addition to the offsets, we need to specify the size
and scale of the grid in $x$ and $y$, as well as the origin of the
grid.

\begin{datatype}
typedef struct {
    int nX;                             ///< Number of elements in x
    int nY;                             ///< Number of elements in y
    double x0;                          ///< Position of 0,0 corner on focal plane 
    double y0;                          ///< Position of 0,0 corner on focal plane
    double xScale;                      ///< Scale of the grid
    double yScale;                      ///< Scale of the grid
    double **x;                         ///< The grid of offsets in x
    double **y;                         ///< The grid of offsets in y
} psFixedPattern;
\end{datatype}

The constructor for \code{psFixedPattern} shall be:
\begin{prototype}
psFixedPattern *psFixedPatternAlloc(double x0,        double y0, 
                                    double xScale,    double yScale,
                                    const psImage *x, const psImage *y);
\end{prototype}
Here, \code{x0}, \code{y0}, \code{xScale} and \code{yScale} have the
same meaning as in the \code{psFixedPattern} structure.  Note that the
values of the fixed pattern offsets are specified as images, the
values from which need to be copied into the \code{double **x} and
\code{double **y} of \code{psFixedPattern}, and that the number of
elements may be derived from the size of the images.

\tbd{Usage of this type is not clear, and awaits prototyping --- do not
worry about coding this in detail yet.}

