%%% $Id: psLibSDRS.tex,v 1.336.2.1 2005-09-19 22:33:44 jhoblitt Exp $
\documentclass[panstarrs,spec]{panstarrs}

% basic document variables
\title{Pan-STARRS PS-1 Image Processing Pipeline Library}
\subtitle{Supplementary Design Requirements}
\shorttitle{PSLib SDRS}
\author{Eugene Magnier, Paul Price, Robert Lupton, Joshua Hoblitt}
\audience{Pan-STARRS PMO}
\group{Pan-STARRS Algorithm Group}
\project{Pan-STARRS Image Processing Pipeline}
\organization{Institute for Astronomy}
\version{16}
\docnumber{PSDC-430-007}

% \setcounter{tocdepth}{5} % lowest level to be included in toc
% \setlength{\topsep}{-2pt}
  
\begin{document}
\maketitle
\sloppy

% make a margin comment:
% \marginpar{note!}

% -- 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 Mar 29 & Draft \\ \hline
00 & 2004 Apr 1  & First version, sent to MHPCC \\ \hline
01 & 2004 May 19 & Extensive modifications, see Appendix B \\ \hline
02 & 2004 Jun 22 & Incorporation of Bugzilla PRs (up to 69) \\ \hline
03 & 2004 Jul 06 & \\ \hline
04 & 2004 Jul 13 & See Appendix B for a change log. \\ \hline
05 & 2004 Aug 16 & draft for start of cycle 3 \\ \hline
06 & 2004 Aug 19 & revision for cycle 3 \\ \hline
07 & 2004 Sep 07 & final for cycle 3 \\ \hline
08 & 2004 Oct 12 & draft for start of cycle 4 \\ \hline
09 & 2004 Nov 15 & final for cycle 4 \\ \hline
10 & 2004 Nov 30 & update for cycle 4 \\ \hline
11 & 2005 Jan 21 & draft for cycle 5 \\ \hline
12 & 2005 Feb 09 & final for cycle 5 \\ \hline
13 & 2005 Mar 30 & draft for cycle 6 \\ \hline
14 & 2005 Apr 27 & final for cycle 6 \\ \hline
15 & 2005 Jun 15 & draft for cycle 7 \\ \hline
16 & 2005 Sep 13 & final for cycle 8 \\ \hline
\RevisionsEnd

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

\DocumentsInternal
PSDC-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
\DocumentsExternal
Posix Standard & Open Group Based Specifications Issue 6, IEEE Std 1003.1, 2003 \\
\DocumentsEnd

\tableofcontents
\pagebreak 
\pagenumbering{arabic}

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

\section{Introduction and policies}

This document describes the Pan-STARRS Image Processing Pipeline (IPP)
Toolkit Library, PSLib.  Any large software project such as the IPP
benefits from the existence of a library of basic software functions
which can be used throughout the project to simplify programming
tasks.  Among the benefits are the ability to reuse code,
simplification of the testing process, streamlining of the code, and
the isolation and encapsulation of concepts which may be subject to
change.  The component functions of such a library should represent
well-defined, concise operations which can be coded with only a modest
number of lines.  PSLib is a library of basic functions required by
the IPP, and it includes many programming concepts which may be useful
for other software projects, especially those which deal with
astronomical data handling tasks.

PSLib consists of a collection of library function calls --- the
Application Programming Interfaces (APIs) --- and the associated data
structures.  The capabilities provided by PSLib are grouped into the
following areas:
%
\begin{itemize}
\item System Utilities
\item Basic Data Collections
\item Data Manipulation
\item Astronomy-Specific Functions.
\end{itemize}
This list is sorted in a hierarchical order: the later entries depend
on the functions and data types of the earlier entries.  

The installed code base for PSLib consists of header files, the binary
library code, \code{libpslib.a} and the shared-library equivalent,
\code{libpslib.so} (or \code{libpslib.dylib} in the case of OS/X).
Assuming these components have been installed into the library and
search path, PSLib may be used within a program by including the line
\code{#include <pslib.h>} into the C code and linking with
\code{-lpslib.}

This document describes the data structures and details the functions
calls. The specified data structures and functions follow the naming
conventions detailed in the IPP Software Requirements Specification
(PSDC-430-005).  In particular, these coding conventions restrict the
namespace used by the library functions by requiring that all globally
visible symbols start with the two letters \code{ps}.  Further
namespace organization is achieved by encouraging functions to be
named in the form psNounVerbPhrase, where Noun is the data type of
principle relevance and VerbPhase describes the operation applied to
that data type.  For example, the function which copies an image (of
type \code{psImage}) is called \code{psImageCopy().}

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

\subsection{External Libraries}

It is anticipated that many of the functions specified in this
document will be implemented through wrapping external libraries:
\begin{itemize}
\item Many of the matrix functions, some of the polynomial and some of
the minimization functions functions should wrap the GNU Scientific
Library (GSL):

\href{www.gnu.org/software/gsl/}{\tt www.gnu.org/software/gsl};

\item The sort functions should wrap the system \code{qsort} call

\item Some of the Fourier transform functions should wrap the Fastest Fourier
Transform in the West library (FFTW):

\href{www.fftw.org}{\tt www.fftw.org}

\item The FITS functions should wrap the CFITSIO library:

\href{heasarc.gsfc.nasa.gov/docs/software/fitsio}{\tt heasarc.gsfc.nasa.gov/docs/software/fitsio}

\end{itemize}

In addition, some of the functions were implemented in a limited
fashion in the Pilot Project.  Also, RHL has provided functional
prototype code for the memory management, tracing and message logging
functions (some of which need to be revised), and some of the data
containers (doubly-linked lists, hashes, arrays).

\subsection{Threads and Re-entrancy}

Due to current developments in CPU architecture, we must assume that
PSLib will be used in a threaded environment.  However, coding the
library to be thread-safe may have implications for the speed of the
library and for the simplicity of use.  We therefore make the
following policies:

\begin{itemize}
\item The memory management and error stack functions, defined below
  (Sections~\ref{sec:memory} \&~\ref{sec:errors}), must be written to
  be thread-safe, since we cannot risk this crucial area being
  unstable.
\item Re-entrant versions of system calls and external library
  functions should be used.  We expect that these cases are
  sufficiently small that we are prepared to err on the side of
  caution. 
\item The practise of using \code{static} variable to achieve high
  efficiency (e.g., so that subsequent calls do not have to repeat a
  large memory allocation) should be kept to an absolute minimum.
  Where it has been justified (i.e., through code profiling), the
  \code{static} variable must be protected by a ``mutex''.
\item Cross-thread synchronization for PSLib's fundamental datatypes
  (\code{psArray}, \code{psList}, etc.) is left to the end-user
  (although some convenience is provided --- see below).
\end{itemize}

As a convenience to the user in achieving thread-safe operation, each
of the data structures classified as a ``collection'' (i.e.,
\code{psList, psHash, psMetadata, psArray, psPixels, psVector,
psBitSet}) and \code{psImage} shall contain a member, \code{void
*lock}, which provides a place for the user to carry around a mutex or
semaphore.  This is provided so that the user doesn't have to pass
around both the structure and a mutex, or wrap PSLib structures in
their own thread-safe structures that contain a mutex.  PSLib is not
responsible for allocating, setting, checking or freeing the
\code{lock} --- these are entirely the responsibility of the user.
PSLib provides only a place to hang it.  Your mileage may vary.

If these policies become a burden on the processing speed, we can
investigate alternative measures, such as defining specifically
re-entrant versions of select speed-critical functions.

\subsection{Angles}

To maintain consistency throughout the library, angles shall be
specified in radians.  In a small number of cases which we expect will
be used heavily (i.e., \code{psSphereOffset} and trigonometric
operations on vectors), the unit may be specified, which provides
convenience to the user.

\subsection{C++ Compatibility}

All PSLib ``public'' header files should comptable with \code{C++}.
This primarily involves using \code{C} pre-processor directives to
enable \code{C++}'s \code{extern "C"} linkage specification when a
header file is being processed as part of a \code{C++} compilation.

An example of wrapping a header file with an \code{extern "C"} block:

\begin{verbatim}
#ifndef FOO_H
#define FOO_H 1

#ifdef __cplusplus
extern "C" {
#endif

...

#ifdef __cplusplus
}
#endif

#endif // FOO_H
\end{verbatim}

Other coding cautions include avoiding trailing commas in
\code{enum}s, protection of the poison \code{pragma}s used for the
memory system (below), and avoiding the qualifier \code{restrict}.

\subsection{Use of {\tt restrict}}

The \code{restrict} type qualifier in C99 indicates to the compiler
that the memory pointed to by a particular pointer is not also pointed
to by some other pointer.  This allows the compiler to optimise the
code, based on the fact that aliasing is not an issue.  However, the
compiler does not check the accuracy of the assumption on which the
optimisation is made, which can lead to data corruption.

Due to the large number of cross-links in psLib (e.g., the metadata
container keeps two pointers to the data), the assumption behind the
use of \code{restrict} generally will not hold.  Correspondingly, use
of \code{restrict} should be kept to a minimum; it should only be
employed where necessary, and where the assumption will hold.
Restricts should also be avoided because they are incompatible with
\code{C++}.

\subsection{Return values}

In some cases, we have defined functions that return a boolean value.
It is intended that (unless otherwise specified), the function returns
\code{true} if there was no error, and returns \code{false} in the
event of an error.

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

\pagebreak 
\section{System Utilities}

\subsection{Configuration}

It is important to be able to access the version of the code in use.
\code{psLibVersion} shall return the current version of PSLib in use,
as determined from CVS tags.

\begin{prototype}
const char *psLibVersion(void);
\end{prototype}

\subsection{Initialization}

Before certain \code{psTime} functions can be employed, the user must
call \code{psTimeInitialize}.

\subsection{Memory Management}
\hlabel{psMemBlock}

\subsubsection{Introduction}

PSLib needs a level of memory management placed between the operating
system (\code{malloc}/\code{free}) and the high level routines (e.g.\
\code{psMetadataAlloc}).  This layer is in addition to the possibility
that specific heavily used data types may need their own
special-purpose memory managers.  However, since we require that all
user-level objects be allocated via associated \code{Alloc/Free}
functions, we will easily be able to implement such functionality
without impacting the facilities described here.

\subsubsection{Rationale}

We wish to insert our own layer of memory management for a number of
reasons:

\begin{itemize}
\item
  We wish to insulate ourselves from the details of the system-provided
  \code{malloc}.  There is no guarantee that the goals of the system
  architect align with those of the PSLib or the IPP.

\item
  We need at least a wrapper layer which handles failed memory
  requests without requiring the application programmer to check
  every attempted allocation.

\item
  We need to provide a mechanism for tracking and fixing memory leaks.
  While it is possible to do this by linking with external libraries
  (e.g.\ \href{gnu.org}{Electric Fence}), it is convenient to do so
  within the framework provided by PSLib.

\item
  Similarly, we wish to provide convenient hooks to detect and diagnose
  memory corruption.

\item
  While debugging complex scientific code, it is very useful to be
  able to trace a specific data structure as it passes through the
  processing pipeline.

\item
  There may be other features that we wish to add in the future (e.g.
  associating a type with every allocation).
\end{itemize}

\subsubsection{Memory Management}
\label{sec:memory}

In the following sections, we specify the API set and define the
appropriate data structures needed by the PSLib memory management
system in order to meet the requirements specified by the desiderata
listed above.

Within the PSLib memory management system, every allocated memory
block which is provided to the user is bounded by two additional
memory segments.  The segment preceeding the user-memory contains data
describing the allocated block, using the \code{psMemBlock} structure.
The first and last elements of this structure are \code{void} pointers
called \code{startblock} and \code{endblock}, which are assigned a
special value, \code{P_PS_MEMMAGIC}.  The segment following the
user-memory block consists of a single \code{void} pointer, and is
also assigned the special value of \code{P_PS_MEMMAGIC}.  The element
\code{userMemorySize} specifies the number of bytes allocated in the
user block and allows the endpost to be found.

In practice, these bounding memory blocks mean that when a user
requests $N$ bytes of memory, the memory management system in fact
allocates \code{N + sizeof(psMemBlock) + sizeof(void)} bytes, starting
at a particular address, \code{ADDR}.  It then fills in the first
\code{sizeof(psMemBlock)} bytes with the data of the \code{psMemBlock}
structure, and the last \code{sizeof(void)} bytes with
\code{P_PS_MEMMAGIC}.  It returns to the user the pointer corresponding
to the address \code{ADDR + sizeof(psMemBlock)}.  If the memory
management system reallocates a block of memory, it must also allocate
the additional space and fill in the boundary values.  If the memory
management system is given a specific pointer for some operation, it
is able to find the corresponding \code{psMemBlock} by simply
subtracting \code{sizeof(psMemBlock)} from the pointer address.

The purpose of the three boundary markers is to catch corruption and
to act as an aid in low-level debugging.  In the first case, memory
over- and under-run errors are likely to overwrite the special values
in either the leading or trailing boundaries.  The typical situation
is one where the coder mis-counts the range and either fills the data
just before the start of the valid memory or just after the end of the
valid memory.  These actions are likely to alter the boundary-post
values, which can be detected by the memory management system.  In the
second case, hexadecimal dumps of large blocks of memory are easier to
examine if the value of \code{P_PS_MEMMAGIC} is chosen to catch the
eye.  A traditional value for \code{P_PS_MEMMAGIC} is \code{0xdeadbeef}
which is also easily recognized in a dump of the memory table.

The structure \code{psMemBlock} specifies additional information
maintained for each block of allocated memory, and is defined as
follows:
%
%\filbreak
\begin{datatype}
typedef struct psMemBlock {
    const void* startblock;             ///< initialised to p_psMEMMAGIC
    struct psMemBlock* previousBlock;   ///< previous block in allocation list
    struct psMemBlock* nextBlock;       ///< next block allocation list
    psFreeFunc freeFunc;                ///< deallocator.  If NULL, use generic deallocation.
    size_t  userMemorySize;             ///< the size of the user-portion of the memory block
    const psMemId id;                   ///< a unique ID for this allocation
    const char* file;                   ///< set from __FILE__ in e.g. p_psAlloc
    const unsigned int lineno;          ///< set from __LINE__ in e.g. p_psAlloc
    pthread_mutex_t refCounterMutex;    ///< mutex to ensure exclusive access to reference counter
    psReferenceCount refCounter;        ///< how many times pointer is referenced
    bool persistent;                    ///< marks if non-user persistent data like error stack, etc.
    const void* endblock;               ///< initialised to p_psMEMMAGIC
} psMemBlock;

typedef void (*psFreeFunc)(void* ptr);
typedef unsigned long psMemId;
typedef unsigned long psReferenceCount;
\end{datatype}
%
The PSLib memory management system must maintain the collection of
allocated memory blocks.  The entries \code{previousBlock} and
\code{nextBlock} point to the previous and next memory blocks
allocated by the memory management system (if they exist, or else
\code{NULL}) and are used to scan through memory blocks as a linked
list.

The element \code{freeFunc} specifies the deallocator associated with
a specific block of memory.  If this element is \code{NULL}, the basic
deallocator (\code{psFree}) is used and the memory block must not be a
rich data structure which requires additional freeing functionality
(otherwise memory leaks will occur).

The element \code{id} in the structure is a sequential memory block
ID.  The memory management system must maintain an internal memory
block ID counter from which a new ID may be supplied to each newly
allocated block of memory and saved in the element \code{id}.  This ID
should also be the key to the memory block in the memory block table.

The two entries \code{file} and \code{lineno} are set to the line
number and file at which the memory was originally allocated.  This is
most easily implemented by the use of the C preprocessor macros
\code{__LINE__} and \code{__FILE__}.  For this reason, we specify
below that the basic memory managment functions be implemented as
preprocessor macros which wrap the intrinsic C level functions.

The element \code{refCounterMutex} is a mutex used to limit access to
the reference counter (below) to a single thread at a time.

The element \code{persistent} indicates if the memory is to be used by
some non-user persistent data.  This allows memory marked as
\code{persistent} to be ignored when identifing memory leaks with
\code{psMemCheckLeaks}.  Private functions shall be provided to test
and set the value of \code{persistent}:
\begin{prototype}
void p_psMemSetPersistent(psPtr ptr, bool value);
bool p_psMemGetPersistent(psPtr ptr);
\end{prototype}


The \code{psMemBlock} structure element \code{refCounter} is provided
so APIs may cleanly manage multiple references to the same block of
memory.  As discussed below, the basic free function, \code{psFree},
is specified to free the memory block only if the reference counter is
set to 1.  See the discussion in section~\ref{sec:free} for an example
of the usage.  Usage of this feature is strongly encouraged, but not
enforced by the memory management system.

In order to trace double frees and other memory errors, the memory
block reference is not automatically deleted when the assocated memory
is deleted.  Rather, the \code{psMemBlock} data and the endpost data
are left behind.  If \code{userMemorySize} is 0, then the memory block
has been freed.  This state must be enforced by \code{psFree}.  This
behavior, in which the associated \code{psMemBlock} is retained, is
only provided if the code is compiled with the preprocessor variable
\code{PS_MEM_DEBUG} defined.

\paragraph{Use of Persistent Memory}

The \code{persistent} member of the \code{psMemBlock} is provided as a
convenience to the end-user of psLib --- it allows him to check for
memory leaks in his own code, without going to the trouble of freeing
memory allocated for psLib core functions.  It also allows the use of
psLib core functions (e.g., \code{psLogMsg, psTrace}) within the
function that is called by \code{psMemCheckLeaks}, without the memory
allocated by those psLib core functions being identified as a leak.

For these reasons, the \code{persistant} feature should only be used
as necessary.  In order to allow proper testing of psLib, all
components that employ any \code{persistent} memory shall provide a
function that frees all of its \code{persistent} memory (and only that
persistent memory that belongs to that component), all of which should
be called from \code{psLibCleanup}.

\subsubsection{APIs for Allocating and Freeing}

PSLib must provide the following APIs to create and destroy memory
blocks:
%
\begin{prototype}
psPtr psAlloc(size_t size);
psPtr psRealloc(psPtr ptr, size_t size);
void psFree(psPtr ptr);
\end{prototype}
%
From the user's perspective, the functions \code{psAlloc},
\code{psRealloc}, and \code{psFree} have identical semantics to the
standard C library functions \code{malloc}, \code{realloc}, and
\code{free}.  In fact, these functions should be implemented as C
preprocessor macros which call the following private functions:
%
\begin{prototype}
psPtr p_psAlloc(size_t size, const char *filename, unsigned int lineno);
psPtr p_psRealloc(psPtr ptr, size_t size, const char *filename, unsigned int lineno);
void p_psFree(psPtr ptr, const char *filename, unsigned int lineno);
\end{prototype}
%
In these function calls, \code{size} is the number of bytes required
to be allocated, \code{ptr} is the pointer to the allocated memory
block, \code{file} is the file containing the calling function (set by
\code{__FILE__}), and \code{line} is the calling function line number
in the source-code file (set by \code{__LINE__}).  Because the user
should only see the preprocessor versions (ie, \code{psAlloc}), we do
not distinguish between these two in the following discussion.

In order to enforce the use of the PSLib versions, the header file
must take steps to ensure that code calling the functions
\code{malloc}, \code{calloc}, \code{realloc}, or \code{free} will not
compile.  This may be achieved by defining preprocessor macros which
mask these functions with invalid statements (\eg{} \code{#define
malloc(S) for}).  In exceptional cases, such as the memory management
system itself, programmers may choose to override this prohibition by
defining the symbol \code{PS_ALLOW_MALLOC}.  Application code will
call \code{p_psAlloc}, \code{p_psRealloc}, or \code{p_psFree} via the
macros defined above.

The above restrictions must apply strictly to the code within PSLib, and other
dependent software, such as e.g., PSModules, may choose to enforce the
restrictions as well.  Software which employs PSLib should be allowed to
optionally employ this enforcement or not, depending on which version of the
global include file is chosen.
  
The functions \code{psAlloc} and \code{psRealloc} must never return a
\code{NULL} pointer. If they are unable to provide the requested
memory they must attempt to obtain the desired memory by calling the
routine registered by \code{psMemExhaustedCallbackSet} (see
\S\ref{secMemAdvanced}), and if still unsuccessful, call
\code{psAbort()}.  The same behavior is true for constructors of rich
structures, with names of the form \code{psFooAlloc}.

Note that we have not specified an equivalent of the \code{calloc}
function.  The \code{calloc} function provides two aspects which
\code{malloc} originally did not: aligned memory and inialization.
Neither of these are required: under POSIX, \code{malloc} is required
to be aligned.  Also, for all structures it is necessary to explicity
define the initialization independently since a byte value of 0 is
usually insufficient.

\begin{prototype}
void psMemSetDeallocator(psPtr ptr, psFreeFunc freeFunc);     
psFreeFunc psMemGetDeallocator(const psPtr ptr);     
\end{prototype}

A free function handles any deallocation procedures of an object
referred to by a pointer, excluding the freeing of the memory block of
the pointer reference itself.  If the pointer refers to an object that
references other memory blocks this free function would insure the
freeing of any encapsulated memory block references. For example,
\code{psList} references a series of node objects of a linked list, so
its free function would free these node objects and the node object's
free function would free the data element references they may contain.

The fucntion \code{psMemSetDeallocator} is used to associated a free
function to a memory block, while \code{psMemGetDeallocator} retrieves
the last free function set.  To remove a free function from a memory
block, the \code{psMemSetDeallocator} should be invoked with
\code{NULL}.  If no free function is set, \code{psMemGetDeallocator}
shall return \code{NULL}.

\paragraph{Negative allocations}

Note that we have specified that the memory size is unsigned
(\code{size_t}), so that we can address the full range of memory that
the architecture will allow, and to match the behaviour of the system
\code{malloc}.  This creates the potential for problems if a negative
value is inadvertently passed to \code{psAlloc} --- it will be
interpreted as a very large positive value.  To guard against this, we
specify that \code{psAlloc} must check that the allocation is less
than \code{PS_MEM_LIMIT} (a preprocessor variable).

For array-like collections (specifically, \code{psArray},
\code{psPixels}, \code{psVector}, and \code{psImage}) we allow the
user to refer to a negative index in the accessor (e.g.,
\code{psArrayGet}) to mean address from the end.  Consequently, the
number of elements in structures should be signed (in order to be able
to access the full range of allocated values).  It is the
responsibility of these structure allocators (e.g.,
\code{psArrayAlloc}) to check that the requested number of elements is
not negative (calling \code{psAbort} otherwise).  All other allocators
shall simply use \code{size_t} where the number of elements is needed
(saving the trouble of checking before passing to psAlloc).

\subsubsection{Callback Routines}

The PSLib memory management system uses callback functions to handle
certain errors and trace conditions.  The callbacks are registered by
the programmer and called by the basic memory management functions if
needed.  The four callbacks currently defined are called in the
following situations:
%
\begin{itemize}
\item when insufficient memory is available (\code{psMemExhaustedCallback})
\item when a memory block is found to be corrupted (\code{psMemProblemCallback})
\item when a specified memory ID is allocated (\code{psMemAllocateCallback})
\item when a specified memory ID is freed (\code{psMemFreeCallback})
\end{itemize}
%
The callback functions are defined in terms of specific callback
types, specified below.  The callbacks are set using functions with
names of the form \code{psCallbackSet}.  In all cases, the \code{Set}
routine takes a pointer to the desired callback function and returns a
pointer to the one that was previously installed. The default
functions for each of these callbacks is listed below:
%
\begin{itemize}
\item \code{psMemExhaustedCallbackDefault}
\item \code{psMemProblemCallbackDefault}
\item \code{psMemAllocCallbackDefault}
\item \code{psMemFreeCallbackDefault}
\end{itemize}
%
and have the behavior of immediately returning \code{NULL}, as if the
callback had been skipped.  If the function pointer passed to the
functions above is \code{NULL}, the callback must be reset to the
default callback function.  The named default callbacks may be used in
within a debugger to catch these conditions by breaking when the
functions are called.  We discuss the use of each of the four
callbacks below.

\subsubsubsection{\tt psMemExhaustedCallback}

If not enough memory is available to satisfy a request by
\code{psAlloc} or \code{psRealloc}, these functions attempt to find an
alternative solution by calling the \code{psMemExhaustedCallback}, a
function which may be set by the programmer in appropriate
circumstances, rather than immediately fail.  The typical use of such
a feature may be when a program needs a large chunk of memory to do an
operation, but the exact size is not critical.  This feature gives the
programmer the opportunity to make a smaller request and try again,
limiting the size of the operating buffer.  This callback has the
following form:
%
\begin{datatype}
typedef psPtr (*psMemExhaustedCallback)(size_t size);
\end{datatype}

\begin{prototype}
psMemExhaustedCallback psMemExhaustedCallbackSet(psMemExhaustedCallback func);
\end{prototype}
%
The callback function is called with the attempted size and is
expected to return a pointer to the allocated memory or \code{NULL}.
Until the callback function is set with
\code{psMemExhaustedCallbackSet}, or if this callback is set to NULL,
\code{psAlloc} immediately calls \code{psAbort}.

\subsubsubsection{\tt psMemProblemCallback}

At various occasions, the memory manager can check the state of the
memory stack.  If any of these checks discover that the memory stack
is corrupted, the \code{psMemProblemCallback} is called.  The callback
has the following form:
%
\begin{datatype}
typedef void (*psMemProblemCallback)(psMemBlock *ptr, const char *filename, unsigned int lineno);
\end{datatype}

\begin{prototype}
psMemProblemCallback psMemProblemCallbackSet(psMemProblemCallback func);
\end{prototype}
%
This callback may be used to report the error and other status
information.  No return value is accepted, and no specific operations
are expected.  The callback is for informational purposes only.  Where
practical and efficient, the memory manager must call the routine
registered using \code{psMemProblemCallbackSet} whenever a corrupted block
of memory is discovered.  For example, doubly-freed blocks can be
detected by checking \code{psMemBlock.refCounter}.

\subsubsubsection{\tt psMemAllocCallback \& psMemFreeCallback}

Two private variables, \code{p_psMemAllocID} and
\code{p_psMemFreeID}, can be used to trace the allocation and freeing
of specific memory blocks.  If the first (\code{p_psMemAllocID}) is
set and a memory block with that ID is allocated, the corresponding
callback is called just before memory is returned to the calling
function.  If the second (\code{p_psMemFreeID}) is set and the memory
block with the ID is about to be freed, the corresponding callback is
called just before the memory block is freed.  These variables are
internal and private to the memory manager and should be set using the
following two functions:
%
\begin{prototype}
psMemId psMemAllocCallbackSetID(psMemId id);
psMemId psMemFreeCallbackSetID(psMemId id);
\end{prototype}
%
The corresponding callback functions have the following form:
%
\begin{datatype}
typedef psMemId (*psMemAllocCallback)(const psMemBlock *ptr);
typedef psMemId (*psMemFreeCallback)(const psMemBlock *ptr);
\end{datatype}

and are set with the following functions:

\begin{prototype}
psMemAllocCallback psMemAllocCallbackSet(psMemAllocCallback func);
psMemFreeCallback psMemFreeCallbackSet(psMemFreeCallback func);
psMemId psMemGetId(void);
\end{prototype}
%
The callback functions are called with a pointer to the corresponding
memory block.  The routines \code{psMemFreeCallbackIDSet} and
\code{psMemAllocCallbackIDSet} accept the desired ID value and return the
old value to the user.  The return values of the handlers installed by
\code{psMemAllocCallbackSet} and \code{psMemFreeCallbackSet} are
used to increment the values of \code{p_psMemAllocID} and
\code{p_psMemFreeID} respectively.  For example, a return value of
\code{0} implies that the value is unchanged; if the value is \code{2}
the callback will be called again when the memory ID counter has
increased by two.  This functionality may be useful to check, for
example, every 100th block allocated.  The function \code{psMemGetId}
returns the next identification number to be assigned to a memory
block.  This function can be used to guide the choice of ID set with
the functions above.

\subsubsection{Memory Tracing and Corruption Checks}
\hlabel{secMemAdvanced}

The PSLib memory management system includes features to facilitate
tracing the memory allocation and freeing process and to debug memory
errors in the calling code.  The types and function prototypes for
this part of the memory API are shown below.
%
\begin{prototype}
int psMemCheckLeaks(psMemId id0, psMemBlock ***array, FILE *fd, bool persistence);
int psMemCheckCorruption(bool abort_on_error);
\end{prototype}
%
The routine \code{psMemCheckLeaks} may be used to check for memory
leaks. The return value is the number of blocks that have been
allocated but not freed.  Only blocks with \code{psMemBlock.id}
greater than \code{id0} are checked; this allows the user to ignore
blocks allocated by initialization routines.

If the argument \code{array} is non-\code{NULL}, then \code{**array}
is set to an array of pointers to \code{psMemBlock} when the function
returns.  These pointers represent the blocks which have been
allocated but not freed.  It is the caller's responsibility to free
this array with \code{psFree}. Also note that \code{**array} should be
\code{NULL} (or not point to allocated memory) upon entering the call
or the corresponding memory reference will be lost.

If the argument \code{fd} is non-\code{NULL}, a one-line summary
of each block that has been allocated but not freed is written to that
file descriptor.

If \code{persistence} is \code{false}, then only those
\code{psMemBlock}s not marked as \code{persistent} shall be considered
as memory leaks.  If \code{persistence} is \code{true}, then all
\code{psMemBlocks} shall be considered as memory leaks.

The routine \code{psMemCheckCorruption} checks the entire heap for
corruption, calling the routine registered with
\code{psMemProblemCallbackSet} for each block detected as being corrupted.
The return value is the number of corrupted blocks detected. If the
argument \code{abort_on_error} is true, \code{psMemCheckCorruption}
must call \code{psAbort} as soon as memory corruption is detected.

\subsubsection{Reference Counting}
\hlabel{secMemRefcounter}

As mentioned above, the memory management system provides the
\code{refCounter} element in \code{psMemBlock} to allow for the
management of multiple references to the same block of memory.
External APIs which make use of this structure must increment the
counter for every additional reference to an allocated memory block,
and decrement it when those references are removed.  The memory
management routines respect the use of the \code{refCounter} field:
\code{psFree} will not free a block for which \code{refCounter} is not
1, but shall decrement the \code{refCounter}, and \code{psAlloc} will
initialize the field to 1.  However, these functions do not (and
cannot practically) enforce the use of the counters; this is a
requirement of external APIs which intend to use the feature.

Several APIs are provided to manage the reference counters.  These
APIs are:
%
\begin{prototype}
psReferenceCount psMemGetRefCounter(const psPtr ptr);
psPtr psMemIncrRefCounter(const psPtr ptr);
psPtr psMemDecrRefCounter(psPtr ptr);
\end{prototype}
%
The functions all take a pointer to the start of a user block of
memory.  The first simply returns the value of the reference counter.
If \code{vptr} is \code{NULL}, this function must return a value of
NULL.  The next two functions increment or decrement the reference
counter, returning the pointer which was passed in. These functions
must validate the memory pointer by determining the corresponding
\code{psMemBlock.id} and checking for consistency in the internal
memory block table (the table pointer for \code{psMemBlock.id} must be
in the valid range and must correspond to the address of the
\code{psMemBlock}).

\subsubsection{Relation of Memory Management to Structures}
\label{sec:free}

Within PSLib and throughout the Pan-STARRS project, we specify a
variety of rich data structures.  The IPP Software Requirements
Specification states that structures should be defined with
corresponding constructors and destructors.  The destructors are
private functions used only by the memory management system.
Instances of, for example, \code{psSomeType} should be constructed using
\code{psSomeTypeAlloc()} calls, and are destroyed using the basic
\code{psFree} function, which calls \code{p_psSomeTypeFree()} to free the
components of the structure, but leaves the task of freeing the
structure to \code{psFree}.  The allocator will allocate the required
memory with \code{psAlloc} and increment the appropriate
\code{refCounter}.

The existence of complicated structures which include pointers to
other structures require that we lay out a rule regarding destructors
and reference counters.  Simply put, \textit{the destructor for every
structure should only free the structure if the \code{refCounter ==
1}; otherwise, it decrements the reference counter and returns.}

\subsubsection{Conventions adopted for pointers}

Only pointers to memory allocated with the PSLib memory functions are
compatible with the various PSLib container types (e.g., \code{psList,
psMetadata}), because the functions working with the container types
search for the attached \code{psMemBlock}.  If a pointer to memory
allocated with another memory system (e.g., the system \code{malloc}),
or generated by offsetting from another pointer that was allocated
with \code{psAlloc}, is used with PSLib, the PSLib functions would
falsely determine that memory is corrupted, because of the missing
\code{psMemBlock}.

To pilot our way through the potential confusion, instead of calling
all pointers (of unspecified type) a ``\code{void*}'', we adopt a
convention, both in this document and to be used in the source, of
referring to a pointer that has a \code{psMemBlock} attached as a
\code{psPtr}:
\begin{datatype}
typedef void* psPtr;
\end{datatype}

For the same reason, we also adopt a convention of referring to a string
that has a \code{psMemBlock} attached as a \code{psString}:
\begin{datatype}
typedef char* psString;
\end{datatype}

That is, \code{psPtr} is used in cases where the function requires the
use of a pointer of unspecified type allocated with \code{psAlloc},
and \code{psString} is used in cases where the function requires the
use of a \code{char*} allocated with \code{psAlloc} (e.g., via
\code{psStringCopy} or \code{psStringNCopy}).

\subsubsection{Strings}

The use of strings within the PSLib memory management system requires
that they be allocated with \code{psAlloc}.  This means that
substrings cannot be placed on containers such as \code{psArray} or
\code{psMetadata}, since these check for the existence of the attached
\code{psMemBlock}.  To get around this, we specify functions that copy
a string into \code{psAlloc}-ed memory:
\begin{prototype}
psString psStringCopy(const char *string);
psString psStringNCopy(const char *string, unsigned int nChar);
\end{prototype}

\code{psStringCopy} shall perform a deep copy of the \code{string}
into \code{psAlloc}-ed memory.  \code{psStringNCopy} shall do the
same, up to a maximum of \code{nChar} characters.

We also specify two useful functions:

\begin{prototype}
ssize_t psStringAppend(char **dest, const char *format, ...)
ssize_t psStringPrepend(char **dest, const char *format, ...)
\end{prototype}

\code{psStringAppend} shall append, according to the \code{format}
values into the \code{dest} string.  \code{dest} shall be allocated if
\code{NULL}, and the length of the new string (excluding the
terminator) shall be returned.  \code{psStringPrepend} shall do
similarly, except it shall prepend to the \code{dest} string.

\subsubsection{Type information}

We need to be able to specify two different sorts of types.

The first, for use with math structures (\S\ref{sec:mathStructures}),
defines a numerical type; e.g., short integer, double-precision
floating point, etc:

\begin{datatype}
typedef enum {
    PS_TYPE_S8,                        ///< Character
    PS_TYPE_S16,                       ///< Short integer
    PS_TYPE_S32,                       ///< Integer
    PS_TYPE_S64,                       ///< Long integer
    PS_TYPE_U8,                        ///< Unsigned character
    PS_TYPE_U16,                       ///< Unsigned short integer
    PS_TYPE_U32,                       ///< Unsigned integer
    PS_TYPE_U64,                       ///< Unsigned long integer
    PS_TYPE_F32,                       ///< Floating point
    PS_TYPE_F64,                       ///< Double-precision floating point
    PS_TYPE_C32,                       ///< Complex numbers consisting of floats
    PS_TYPE_C64,                       ///< Complex numbers consisting of doubles
    PS_TYPE_BOOL                       ///< Boolean value
} psElemType;
\end{datatype}


The second, primarily for use with the metadata
(\S\ref{sec:metadata}), defines a data structure; e.g., list, array,
FITS file:

\begin{datatype}
typedef enum {                         ///< type of item.data is:
    PS_DATA_S32  = PS_TYPE_S32,        ///< psS32
    PS_DATA_F32  = PS_TYPE_F32,        ///< psF32
    PS_DATA_F64  = PS_TYPE_F64,        ///< psF64
    PS_DATA_BOOL = PS_TYPE_BOOL,       ///< psBool
    PS_DATA_STRING = 0x10000,          ///< String (char *)
    PS_DATA_ARRAY,                     ///< psArray
    PS_DATA_BITSET,                    ///< psBitSet
    PS_DATA_CUBE,                      ///< psCube
    PS_DATA_FITS,                      ///< psFits
    PS_DATA_HASH,                      ///< psHash
    PS_DATA_HISTOGRAM,                 ///< psHistogram
    PS_DATA_IMAGE,                     ///< psImage
    PS_DATA_KERNEL,                    ///< psKernel
    PS_DATA_LIST,                      ///< psList
    PS_DATA_LOOKUPTABLE,               ///< psLookupTable
    PS_DATA_METADATA,                  ///< psMetadata
    PS_DATA_METADATAITEM,              ///< psMetadataItem
    PS_DATA_MINIMIZATION,              ///< psMinimization
    PS_DATA_PIXELS,                    ///< psPixels
    PS_DATA_PLANE,                     ///< psPlane
    PS_DATA_PLANEDISTORT,              ///< psPlaneDistort
    PS_DATA_PLANETRANSFORM,            ///< psPlaneTransform
    PS_DATA_POLYNOMIAL1D,              ///< psPolynomial1D
    PS_DATA_POLYNOMIAL2D,              ///< psPolynomial2D
    PS_DATA_POLYNOMIAL3D,              ///< psPolynomial3D
    PS_DATA_POLYNOMIAL4D,              ///< psPolynomial4D
    PS_DATA_PROJECTION,                ///< psProjection
    PS_DATA_SCALAR,                    ///< psScalar
    PS_DATA_SPHERE,                    ///< psSphere
    PS_DATA_SPHEREROT,                 ///< psSphereRot
    PS_DATA_SPLINE1D,                  ///< psSpline1D
    PS_DATA_STATS,                     ///< psStats
    PS_DATA_TIME,                      ///< psTime
    PS_DATA_VECTOR,                    ///< psVector
    PS_DATA_UNKNOWN,                   ///< Other data of an unknown type
    PS_DATA_METADATA_MULTI             ///< Used internally for metadata; not a 'real' type
} psDataType;
\end{datatype}

Here we have not included every type of structure used in PSLib, but
only those we expect will be frequently carried around in containers.
We include \code{PS_DATA_S32, PS_DATA_F32, PS_DATA_F64} and
\code{PS_DATA_BOOL} since the metadata use these, and we set them to
the same values as in \code{psElemType}.

The other values are offset from these ``primitive'' types so that
they can be identified easily.  A macro,
\code{PS_DATA_IS_PRIMITIVE(type)}, shall return \code{true} if the
type is one of the numerical data types (S32, F32, F64, bool).  In
such a case, the data value is directly available from the metadata.
Otherwise, a pointer to the data is available.

\code{PS_DATA_METADATA_MULTI} is used by the metadata, so the user
should not use this type.

\subsubsection{Type checking}

Several of the collections contain data using a \code{void*} pointer.
This makes it difficult to ensure that data coming off a collection is
of a particular, desired type.  Nevertheless, there is a way of
identifying the type of a pointer pulled off a collection --- the
registered deallocator function stored in the \code{psMemBlock}.  We
specify functions below, one for each of the major types in psLib,
that check the type on a particular pointer, returning \code{true} if
the \code{ptr} matches the desired type, as determined from the
registered deallocator function.  These functions may be implemented
as macros or inline functions if that is deemed convenient.

\begin{prototype}
bool psMemCheckArray(psPtr ptr);
bool psMemCheckBitSet(psPtr ptr);
bool psMemCheckCube(psPtr ptr);
bool psMemCheckFits(psPtr ptr);
bool psMemCheckHash(psPtr ptr);
bool psMemCheckHistogram(psPtr ptr);
bool psMemCheckImage(psPtr ptr);
bool psMemCheckKernel(psPtr ptr);
bool psMemCheckList(psPtr ptr);
bool psMemCheckLookupTable(psPtr ptr);
bool psMemCheckMetadata(psPtr ptr);
bool psMemCheckMetadataItem(psPtr ptr);
bool psMemCheckMinimization(psPtr ptr);
bool psMemCheckPixels(psPtr ptr);
bool psMemCheckPlane(psPtr ptr);
bool psMemCheckPlaneDistort(psPtr ptr);
bool psMemCheckPlaneTransform(psPtr ptr);
bool psMemCheckPolynomial1D(psPtr ptr);
bool psMemCheckPolynomial2D(psPtr ptr);
bool psMemCheckPolynomial3D(psPtr ptr);
bool psMemCheckPolynomial4D(psPtr ptr);
bool psMemCheckProjection(psPtr ptr);
bool psMemCheckScalar(psPtr ptr);
bool psMemCheckSphere(psPtr ptr);
bool psMemCheckSphereRot(psPtr ptr);
bool psMemCheckSpline1D(psPtr ptr);
bool psMemCheckStats(psPtr ptr);
bool psMemCheckTime(psPtr ptr);
bool psMemCheckVector(psPtr ptr);
\end{prototype}

For user convenience, we also specify a one-stop shop which simply
executes the appropriate \code{psMemCheckWhatever} function according
to the \code{type}:

\begin{prototype}
bool psMemCheckType(psDataType type, psPtr ptr);
\end{prototype}

The reason for having functions that check the type instead of
returning the type is that the check is quick and the result is clean,
while returning the type involves descending through a case statement.

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

%%%
%%% Not sure we need these --- all we do is wrap pthread_mutex, so why bother?
%%%

%% \subsection{Locks}
%% 
%% We define the following conveniences:
%% \begin{datatype}
%% typedef pthread_mutex_t psMutex;
%% \end{datatype}
%% 
%% \begin{prototype}
%% psMutex *psMutexAlloc(void);
%% bool psMutexLock(psMutex *mutex);
%% bool psMutexUnlock(psMutex *mutex);
%% \end{prototype}
%% 
%% \code{psMutex} simply wraps a POSIX thread mutex (this is necessary in
%% order to use the PS memory management system).
%% 
%% \code{psMutexAlloc} shall allocate memory for a \code{psMutex}, and
%% initialise the mutex.  \code{psMutexLock} shall lock the \code{mutex},
%% and \code{psMutexUnlock} shall unlock the \code{mutex}.
%% 
%% These functions, in the interests of speed, should be implemented as
%% preprocessor macros.
%% 
%% These functions and the \code{void *lock} in the collection structures
%% and  \code{psImage} should  only  be defined  if \code{_REENTRANT}  is
%% defined; otherwise the functions should be poisoned.

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

\subsection{Tracing and Logging}

This section defines the \PS{} Tracing and Logging APIs; the former
refers to debug information that we wish to be able to turn on and off
without recompiling (it will \emph{not} be available in production
code); the latter means information about the processing that must be
collected and saved, even in the production system.  We envision that
extensive use will be made of \code{psTrace} throughout the \PS{}
code.

\subsubsection{Tracing APIs}
\hlabel{psTrace}

A code-tracing facility should allow the programmer to place messages
in the code which, when called, will print some useful information
about the containing block of code.  Ideally, different messages may
be specified to have different levels (of severity or interest).  For
a given run of the program, the level of interest should be set to
provide more or less feedback, depending on the needs of the
programmer.  In a typical situation, low-level messages would be
placed generously throughout the code, indicating the flow of the
program.  Higher-level messages would be placed in a limited number of
special locations, such as the start of major code segments or where a
particularly unusual condition is met.  Top-level messages would be
placed in code triggered under serious error conditions.  A normal run
of the program would have the trace messages printed only for the
top-level.  If the user needs to dig deeper into the code, the trace
level should be set lower, and the more detailed messages could be
examined.  In a case of a serious, poorly-understood problem with the
code, the trace threshold would be placed to the bottom and the
lowest-level step-by-step messages would be printed.

The PSLib tracing facility provides the above functionality, along
with the ability to assign different trace levels to messages from
different software components.  Each trace message, when placed in the
code, is assigned to be part of a specific tracing 'facility', defined
in more detail below.  The trace level for that specific message is
also set when the message is placed.  Each facility may have its trace
level set independently.  Thus, it is possible to request detailed
trace output for one facility while minimizing the verbosity of the
trace output from the rest of the program.

The trace facilities consist of a hierarchy of names.  A trace
facility is defined by a string consisting of words separated by dots,
with a hierarchy equivalent to that of UNIX directory names.  The
top-level facility is simply \code{'.'} (one dot).  The next level
would be \code{'.A'}, followed by \code{.A.B}, and so on.  The
relationship is seen in two ways.  First, a facility inherits the
trace level of its parent unless explicitly specified.  Second, the
hierarchy is used to format the listing of the trace facilities so
that they are easy to read.  The first of these rules provides a
mechanism to define the default trace levels for any facility even if
it has not been registered explicitly since all named facilities are
implicitly children of the top level facility (\code{.}).  The second
rule is simply an organizational technique to make the listing of
facility information clear.  In specifying a facility, the leading
dot shall be optional, as a convenience to the user.

The API to place a trace message in the code, and simultaneously set
its trace level and facility, is:
%
\begin{prototype}
void psTrace(const char *facil, int level, const char *format,...);
void psTraceV(const char *facil, int level, const char *format, va_list ap);
\end{prototype}
% 
where the \code{format} argument is a printf-style formatting code
followed by possible arguments to that formatting statement, to be
implemented using the \code{vprintf} functions.  This command
specifies the name of the facility to which the message belongs
(\code{facil}), the trace level for this message in that facility
(\code{level}) and the message itself.  The \code{psTraceV} version of
the command accepts a \code{va_list} argument list rather than a
variable number of arguments.

The trace level for any facility may be set at any time with the
following function:
%
\begin{prototype}
int psTraceSetLevel(const char *facil, int level);
\end{prototype}
% 
where \code{level} specifies the current trace level for the facility
named by \code{facil}.  The function returns the previous trace level
for that facility.  The currently defined trace level for a given
facility may be determined by the function:
%
\begin{prototype}
int psTraceGetLevel(const char *facil);
\end{prototype}
% 
which returns the trace level of the named facility following the
rules specified above.  A specified trace message (identified by
\code{psTrace}) must be printed if and only if
\code{psTraceGetLevel(facil)} returns a value greater than or equal to
the value of \code{level} for that message.  That is, a larger
number for the trace level corresponds to lower-level statements, and
hence is more verbose.

PSLib includes a utility function for examining the current tracing
levels of all facilities: 
%
\begin{prototype}
void psTracePrintLevels(void);
\end{prototype}
%
This function prints the hierarchy of trace facilities along with the
current trace level for each facility.  For example, a particular
program may have a few facilities defined, along with their trace
levels.  A call to \code{psTracePrintLevels} may produce a listing
which appears as:
\begin{verbatim}
.                        0
 .IPP                    0
 .IPP.debias             1
 .IPP.flatten            1
  .IPP.flatten.divide    2
  .IPP.object.findpeak   1
  .IPP.object.getsky     1
\end{verbatim}

Considering the same program, the programmer might place a variety of
trace messages throughout the \code{flatten} portion of the code with
different types of messages, such as:
%
\begin{verbatim}
psTrace("IPP.flatten", 2, "starting flatten function\n");
psTrace("IPP.flatten", 0, "flat-field image is \%s\n", filename);
psTrace("IPP.flatten.divide", 2, "doing the divide\n");
psTrace("IPP.flatten.divide", 3, "trying the loop, i = \%d\n", i);
psTrace("IPP.flatten.divide", 1, "got an invalid pixel value (NaN) at \%d,\%d\n");
psTrace("IPP.flatten.divide", 2, "divide is done\n");
\end{verbatim}
%
Under the trace levels set above, if the code actually reached each of
these trace messages, the following messages would be printed:
%
\begin{verbatim}
flat-field image is foo.fits
  doing the divide
  got an invalid pixel value (NaN) at 500,20
  divide is done
\end{verbatim}
%
while these two would not be printed because their facility level was
too low:
%
\begin{verbatim}
  starting flatten function
   trying the loop, i = 0   
   trying the loop, i = 1   
   trying the loop, i = 2   
...
\end{verbatim}
%

The availability of the tracing facility at run-time, must be decided
at compilation: If the C pre-processor macro \code{PS_NO_TRACE} is
defined, all trace code must be replaced by empty space so that none
of the code is compiled.  This can be implemented via macro front-ends
to private versions of the user APIs.  In addition, a function
\code{bool psTraceReset(void)} will free memory used by \code{psTrace}
functions, effectively resetting all trace levels to 0.

The trace may optionally be written to a file or other output
destination with \code{psTraceSetDestination}:
\begin{prototype}
void psTraceSetDestination(int fd);
\end{prototype}
If the \code{fd} is zero, then tracing is disabled.  Otherwise, the
trace is sent to the specified file descriptor.  As a convenience, the
following are defined:
\begin{datatype}
enum {
    PS_TRACE_TO_NONE = 0,                ///< turn off all traces
    PS_TRACE_TO_STDOUT = 1,              ///< trace to system's stdout
    PS_TRACE_TO_STDERR = 2,              ///< trace to system's stderr
};
\end{datatype}
%
This arrangement mirros the file descriptors for standard input,
output and error.  A call to \code{psTraceSetDestination}
automatically closes an open file descriptor.

The corresponding function
\begin{prototype}
int psTraceGetDestination(void);
\end{prototype}
returns the current trace destination file descriptor.  If the
destination has not been defined by the user, the descriptor for
\code{stdout} is returned.

The trace output format is controlled with the function:
%
\begin{prototype}
bool psTraceSetFormat(const char *format);
\end{prototype}
%
which expects a string consisting of the letters \code{H} (host),
\code{L} (level), \code{M} (message), \code{N} (name), \code{F}
(file:name), and \code{T} (time).  The default is \code{THLNM}, which
produces trace messages in the form:
\begin{verbatim}
YYYY-MM-DD hh:mm:ssZ | hostname | L | name | file:line
    The message goes here
    and is indented by 4 spaces.
\end{verbatim}
where \code{YYYY}, \code{MM}, \code{DD}, \code{hh}, \code{mm}, and
\code{ss} are the year, month (Jan is 01), day of the month, hours
(0--23), minutes, and seconds when the trace message was received.
Note that the timestamp is in ISO order, and that the timezone is GMT
(hence the \code{Z}).  The \code{hostname} is returned by
\code{gethostname}, \code{L} is the numerical level. The other two
fields, \code{name} and \code{msg}, are the facility name and the
complete message provided to \code{psTrace}.  The \code{msg} is placed
on a new line (allowing the \code{name} to fill the rest of the
previous line), with each line indented by 4 spaces.  An example
message is:
%
\begin{verbatim}
2004:02:24 20:14:18Z | alibaba.IfA.Hawaii.Edu | I | example.utils.helloWorld | example.c:20
    Hello world,
    it's me calling.
\end{verbatim}
%
The possible order of the format entries is fixed and not determined
by the order of the letters used in \code{psTraceSetFormat}.  Selecting
an output format with fewer than the complete set of 5 entries simply
removes those entries from the output messages.

Specifying a \code{format} of \code{NULL} turns off logging (equivalent
to calling \code{psTraceSetDestination(PS_TRACE_TO_NONE)}, whereas if the
\code{format} is \code{""}, then the format reverts to the default.

\subsubsection{Message Logging}
\hlabel{psLogMsg}

Message logging is similar in some respects to tracing.  Like trace
messages, log messages are placed in the code at various locations to
provide output describing the current state of the program.  Like
the PSLib trace facility, a good log facility should have the
capability of associating each message with an importance or severity
level, and at any point, the level for which messages are actually
printed should be set in a flexible manner.   Unlike trace messages,
however, log messages are always part of the code and are available in
the production version as well as in test versions.  

The PSLib logging facility does not include the extensive facility
levels which are provided by the trace facility.  Less control over
the granularity is needed for the log messages than for the trace
messages.  

A log message is placed in the code with the command:
%
\begin{prototype}
void psLogMsg(const char *name, int level, const char *format, ...); 
void psLogMsgV(const char *name, int level, const char *format, va_list ap); 
\end{prototype}
where \code{name} is a word to describe the source of the message,
\code{level} is the severity level of this message, and \code{format}
is a printf-style formatting statement defining the actual message,
and is followed by the arguments to the formatting code.  The second
form, \code{psLogMsgV} is an equivalent command which takes a
\code{va_list} argument.

A log message may have any level specified in the range 0-9, though
the first 4 levels are associated with symbolic names:
%
\begin{datatype}
enum { PS_LOG_ABORT = 0, PS_LOG_ERROR, PS_LOG_WARN, PS_LOG_INFO };
\end{datatype}
%

At any time, the program may set the current log level, the level
above which log messages are ignored, using the function:
%
\begin{prototype}
int psLogSetLevel(int level);           
\end{prototype}
%
This function returns the previous log level.  A specific message
invoked with \code{psLogMsg} is only printed if its value of
\code{level} is less than the current value set by
\code{psLogSetLevel}.  The default log level is set to
\code{PS_LOG_INFO}.

\begin{prototype}
int psLogGetLevel();           
\end{prototype}
%
This function returns the current log level.  

Log messages are sent to the destination most recently set using:
%
\begin{prototype}
bool psLogSetDestination(int fd);      
\end{prototype}
%
If the \code{fd} is zero, then logging is disabled.  Otherwise, the
log is sent to the specified file descriptor.  As a convenience, the
following are defined:
\begin{datatype}
enum {
    PS_LOG_TO_NONE = 0,                  ///< turn off logging
    PS_LOG_TO_STDOUT = 1,                ///< log to system's stdout
    PS_LOG_TO_STDERR = 2,                ///< log to system's stderr
};
\end{datatype}
%
This arrangement mirros the file descriptors for standard input,
output and error.  A call to \code{psLogSetDestination} automatically
closes an open file descriptor.

The corresponding function
\begin{prototype}
int psLogGetDestination();
\end{prototype}
returns the current log destination file descriptor.  If the
destination has not been defined by the user, the descriptor for
\code{stdout} is returned.

The output format is controlled with the function:
%
\begin{prototype}
bool psLogSetFormat(const char *format);
\end{prototype}
%
which expects a string consisting of the letters \code{H} (host),
\code{L} (level), \code{M} (message), \code{N} (name), \code{F}
(file:line), and \code{T} (time).  The default is \code{THLNM}, which
produces log messages in the form:
\begin{verbatim}
YYYY-MM-DD hh:mm:ssZ | hostname | L | name | file:line
    The message goes here
    and is indented by 4 spaces.
\end{verbatim}
where \code{YYYY}, \code{MM}, \code{DD}, \code{hh}, \code{mm}, and
\code{ss} are the year, month (Jan is 01), day of the month, hours
(0--23), minutes, and seconds when the log message was received.  Note
that the timestamp is in ISO order, and that the timezone is GMT
(hence the \code{Z}).  The \code{hostname} is returned by
\code{gethostname}, \code{L} is a character associated with the level
(\code{A}, \code{E}, \code{W}, and \code{I} for \code{PS_LOG_ABORT},
\code{PS_LOG_ERROR}, \code{PS_LOG_WARN}, and \code{PS_LOG_INFO}
respectively. Other levels are represented numerically (\code{5}
etc.). The other two fields, \code{name} and \code{msg}, are the
arguments to \code{psLogMsg}.  The \code{msg} is placed on a new line
(allowing the \code{name} to fill the rest of the previous line),
with each line indented by 4 spaces.  An example message is:
%
\begin{verbatim}
2004:02:24 20:14:18Z | alibaba.IfA.Hawaii.Edu | I | example.utils.helloWorld
    Hello world,
    it's me calling.
\end{verbatim}
%
The possible order of the format entries is fixed and not determined
by the order of the letters used in \code{psLogSetFormat}.  Selecting
an output format with fewer than the complete set of 5 entries simply
removes those entries from the output messages.

Specifying a \code{format} of \code{NULL} turns off logging (equivalent
to calling \code{psLogSetDestination(PS_LOG_TO_NONE)}, whereas if the
\code{format} is \code{""}, then the format reverts to the default.

The following utility opens an output file descriptor for use by the
trace and log facilities. 
\begin{prototype}
int psMessageDestination (const char *dest);
\end{prototype}
%
The destination string consists of a URL in the form
\code{protocol:location}.  The \code{protocol} value may be
\code{file}, to send the log to a local file named by the value of
\code{location}.  Future expansion may allow the logger to send
messages to an IP logger, with a protocol to be defined later.  Three
other special values are allowed for the \code{dest} parameter
(without specifying a protocol): \code{stderr} and \code{stdout},
which return the file descriptors for \code{stderr} and \code{stdout}
respectively, and \code{none} which returns the special descriptor to
turn off logging.

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

\subsection{Error Handling}
\label{sec:errors}
\hlabel{errorStack}

\PS{} errors shall be raised using a function, \code{psError}, with the
caller supplying a component name and error message.  It is desireable
to be able to trace an error through the program so that the events
that led to the error may be quickly and clearly identified.
\code{psError} prints an error message and doesn't abort, but instead
returns the error code.
\begin{prototype}
psErrorCode p_psError(const char *filename, 
                      unsigned int lineno, 
                      const char *func, 
                      psErrorCode code,
                      bool new, 
                      const char *format, ...);
\end{prototype}
\begin{datatype}
#define psError(code, new, format, ...) \
    p_psError(__FILE__, __LINE__, __func__, code, new, format, __VA_ARGS__)
\end{datatype}

\code{psError} is a macro definition that allows the filename, line
number and function name to be inputted to a private function,
\code{p_psError}.  The \code{code} is an enumerated type which lists
the possible \textit{classes} of errors (e.g. \code{PS_ERR_IO}) that
\PS{} code can generate (see section~\ref{psErrorCodes}). The
\code{new} argument takes a boolean which, if \code{true} specifies
that the error was set initially at this location, and if \code{false}
specifies that an error was passed to this location.  Raising new
error should clear the error stack.  The final required argument,
\code{format}, is a \code{printf}-style format that is passed to
\code{psLogMsgV} with code \code{PS_LOG_ERROR}, and component equal to
the concatenation of the file name and the line number, separated by a
colon.  The result of a call to \code{psError} shall be to push an
error onto a stack; this stack is cleared if \code{new} is true, or by
a call to \code{psErrorClear}.

The errors on the error stack are defined as the following:
\begin{datatype}
typedef struct {
    char *name;                         ///< category of code that caused the error
    psErrorCode code;                   ///< class of error (equivalent to errno)
    char *msg;                          ///< the message associated with the error
} psErr;
\end{datatype}

The last error reported is available from \code{psErrorLast}; if no
errors are current, a non-\code{NULL} \code{psErr} shall be returned
with code \code{PS_ERR_NONE}.  Previous errors on the stack shall be
returned by \code{psErrorGet} (a value of \code{0} passed to
\code{psErrorGet} is equivalent to a call to \code{psErrorLast}).
The error stack may be completely cleared with \code{psErrorClear}.
%
\begin{prototype}
unsigned int psErrorGetStackSize(void);
const psErr *psErrorGet(int which);
const psErr *psErrorLast(void);
void psErrorClear(void);
\end{prototype}

\code{psErrorGetStackSize} shall return the number of errors on the
stack.  The entire error stack may be printed to an open file
descriptor by calling \code{psErrorStackPrint} (or
\code{psErrorStackPrintV}); if and only if there are current errors,
the printf-style string \code{format} is first printed to the file
descriptor \code{fd}. In this printout, error codes shall be replaced
by their string equivalents as defined in the next section.  Note that
these are also available in the \code{psErr} structure. The successive
lines of the traceback should be indented by an additional space.
\code{psErrorStackPrintV} must not invoke \code{va_end}.
%
\begin{prototype}
void psErrorStackPrint(FILE *fd, const char *format, ...);
void psErrorStackPrintV(FILE *fd, const char *format, va_list va);
\end{prototype}

Any error \code{code}s less then or equal to \code{PS_ERR_BASE} (see
next section) must be taken to be valid values of \code{errno}, and
\code{psErrorStackPrint} must print the value returned by
\code{strerror} if such error codes are encountered.

The routine \code{psErrorCodeString} returns the string associated
with an error code:
\begin{prototype}
const char *psErrorCodeString(psErrorCode code);
\end{prototype}

\subsubsection{Error Codes}
\hlabel{psErrorCodes}

Both error codes for PSLib and error codes for projects that use PSLib
may be registered.  In the former case, the error codes must be
registered on initialisation (see \code{psLibInit}), whereas in the
latter case, they must be explicitly registered by the programmer.

\paragraph{Registering error codes}

PSLib and any project needed to use PSLib must define the necessary
error codes and associated message strings.  An array of error codes
may be registered with the PSLib error handler using the function:
\begin{prototype}
void psErrorRegister(const psErrorDescription *errors, int errorCode);
\end{prototype}
where the errors are represented internally as follows:
\begin{datatype}
typedef struct {
    psErrorCode code;                  ///< An error code
    const char *description;           ///< the associated description
} psErrorDescription;
\end{datatype}
PSLib internal errors must be registered with the function
\code{psErrorRegister}, which should be called as part of the
program initialization to set up the error codes.  It is left to the
external project to produce their own error registration functions
which must also be called during initialization. There is a clear need
to coordinate the choice of error numbers.  It is expected that error
code ranges for different projects must be managed by the Project
Office within Pan-STARRS.

\paragraph{Error Codes for PSLib}

For ease of maintenance, error codes for PSLib must be defined by an
auxiliary file, conventionally named \file{psErrorCodes.dat}.  The
format of this file must consist of a number of lines, each of the
form:
\begin{verbatim}
NAME [ = value ][,] STRING
\end{verbatim}
where \code{[ = value]} and the comma are optional, and no spaces are
significant except in the STRING.  Comments extend from \code{#} to
the end of the line (except that a \code{\#} must be replaced by
\code{#} and not taken to start a comment). For example,
\begin{verbatim}
#
# This file is used to generate psErrorClasses.h
#
NONE = 0,               not an error; must be 0
BASE = 256,             first value we use; should avoid errno conflicts
UNKNOWN,                unknown error
# This is a comment, and is ignored.
IO,                     I/O error
BADFREE,                bad argument to psFree()
MEMORY_CORRUPTION,      memory corruption detected
\end{verbatim}
The values \code{NONE = 0} and {UNKNOWN} must be present.

A script, called from the Makefiles, must generate two files,
\file{psErrorCodes.h} and \file{psErrorCodes.c} from the input file
\file{psErrorCodes.dat}.  \file{psErrorCodes.h} must define an
enumerated type \code{psErrorCode} with elements \code{PS_ERR_NAME}
and values as specified in \file{psErrorCodes.dat}, e.g.
\begin{datatype}
#if !defined(PS_ERROR_CODES_H)
#define PS_ERROR_CODES_H

typedef enum {
    PS_ERR_NONE = 0,
    PS_ERR_BASE = 256,
    PS_ERR_UNKNOWN,
    PS_ERR_IO,
    PS_ERR_BADFREE,
    PS_ERR_MEMORY_CORRUPTION,
    PS_ERR_N_ERR_CLASSES,
} psErrorCode;
#endif
\end{datatype}

\file{psErrorCodes.c} must define the necessary functions to register
the error codes.

This script will likely be of use to the user, and so it shall be
installed as part of PSLib.


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

\subsection{Abort}

\code{psAbort}, must call \code{psLogMsgV} with a level of
\code{PS_LOG_ABORT}, and then call \code{abort}.

\begin{prototype}
void psAbort(const char *name, const char *format, ...);
\end{prototype}

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

\subsection{Command-line arguments}

The following functions are provided to aid parsing of command-line
arguments:

\begin{prototype}
int psArgumentGet(int argc, char **argv, const char *arg);
bool psArgumentRemove(int argnum, int *argc, char **argv);
\end{prototype}

\code{psArgumentGet} shall return the index of the element in the
argument list, \code{argv} with number of entries \code{argc}, that
matches the specified argument, \code{arg}, or zero if there is no
match.

\code{psArgumentRemove} shall remove from the argument list
(\code{argv} with number of entries \code{argc}) the argument whose
index is \code{argnum}.  The number of entries in the argument list
shall be decremented.  The function shall return \code{true} if the
\code{argnum} is in the argument list, and \code{false} otherwise.

\begin{prototype}
int psArgumentVerbosity(int *argc, char **argv);
\end{prototype}

\code{psArgumentVebosity} shall implement the various verbosity
controls under the following guidelines:
\begin{itemize}
\item \code{-v} switch shall set the log level to 3.
\item \code{-vv} switch shall set the log level to 4.
\item \code{-vvv} switch shall set the log level to 5.
\item \code{-logfmt someFormat} switch shall set the log format to
  \code{someFormat}.
\item \code{-trace facility level} switch shall set the trace level
  for the specified \code{facility} to the specified \code{level}.
\item \code{-trace-levels} switch shall print the trace levels as
  currently set.
\end{itemize}

The above arguments shall be removed from the argument list as they
are processed.  The function shall return the resultant logging level.

\begin{prototype}
bool psArgumentParse(psMetadata *arguments, int *argc, char **argv);
\end{prototype}

\code{psArgumentParse} shall parse the command line arguments
(supplied via \code{argc, argv}) into a metadata container of
\code{arguments}.  The input \code{arguments} shall contain the list
of possible arguments as the keywords providing the default values (of
the appropriate type).  As matching arguments are found on the command
line, the values shall be read into the \code{arguments} metadata,
with the appropriate type.  Multiple values (i.e., when an argument
takes two or more values) are specified using a
\code{PS_DATA_METADATA_MULTI} list (i.e., by specifying a type with
the \code{PS_META_DUPLICATE_OK} flag).  A boolean argument shall be
set to \code{true} by the presence of the argument itself (no value is
specified).  The arguments and their values shall be removed from the
list of command line arguments as they are processed.  The function
shall return \code{false} without modifying the \code{arguments} if
any argument (i.e., a string beginning with a dash; but not a string
that does not start with a dash, which is likely a mandatory argument
for the program; the purpose of this requirement is so that the
default values are preserved) was encountered that is not present in
the \code{arguments}.

\begin{prototype}
void psArgumentHelp(psMetadata *arguments);
\end{prototype}

\code{psArgumentHelp} shall print to \code{stdout} a guide to the
command-line \code{arguments} (containing the same information as for
\code{psArgumentParse}) of the form:
\begin{verbatim}
Optional arguments, with default values:
    -names     (FRED)     This is the comment from the first value
               (JIM)      This is the comment from the second value
               (BOB)      This is the comment from the third value
    -answer    (42)       This is a comment
    -truth     (FALSE)    This is another comment
\end{verbatim}
Here the \code{arguments} contains three keywords, \code{names, magic,
truth}.  \code{names} has three values, \code{FRED, JIM, BOB}, stored
as a \code{PS_DATA_METADATA_MULTI}, each with the comment as shown.

It will be the responsibility of the user to supply information for
the mandatory arguments.  An example follows:
\begin{verbatim}
int main(int argc, char *argv[])
{
    // Parse optional command-line arguments
    psMetadata *arguments = psMetadataAlloc(); // The arguments, with default values
    psMetadataAdd(arguments, PS_LIST_TAIL, "-string", PS_META_STR, "Test string", "SomeString");
    psMetadataAdd(arguments, PS_LIST_TAIL, "-bool", PS_META_BOOL, "Test bool", false);
    psMetadataAdd(arguments, PS_LIST_TAIL, "-int", PS_META_S32 | PS_META_DUPLICATE_OK, "Test integer 1", 1);
    psMetadataAdd(arguments, PS_LIST_TAIL, "-int", PS_META_S32 | PS_META_DUPLICATE_OK, "Test integer 2", 2);
    psMetadataAdd(arguments, PS_LIST_TAIL, "-int", PS_META_S32 | PS_META_DUPLICATE_OK, "Test integer 3", 3);
    psMetadataAdd(arguments, PS_LIST_TAIL, "-float", PS_META_F32, "Test float", 1.234567);
    if (! psArgumentParse(arguments, &argc, argv) || argc != 3) {
	printf("\nName of the program here\n\n");
	printf("Usage: %s INPUT OUTPUT\n\n", argv[0]);
	psArgumentHelp(arguments);
	psFree(arguments);
	exit(EXIT_FAILURE);
    }
    const char *inputName = argv[1];	// Name of input
    const char *outputName = argv[2];	// Name of output
    printf("Success: %s %s\n", inputName, outputName);

    psString string = psMetadataLookupString(NULL, arguments, "-string");
    float floating = psMetadataLookupF32(NULL, arguments, "-float");
    bool boolean = psMetadataLookupBool(NULL, arguments, "-truth");
    printf("String: %s\n", string);
    printf("Float: %f\n", floating);
    printf("Boolean: %d\n", boolean);

    psMetadataItem *intItem = psMetadataLookup(arguments, "-int");
    if (intItem->type == PS_META_MULTI) {
	psList *intMulti = intItem->data.V;
	psListIterator *intIter = psListIteratorAlloc(intMulti, PS_LIST_HEAD, false);
	psMetadataItem *intMultiItem = NULL;
	while (intMultiItem = psListGetAndIncrement(intIter)) {
	    printf("Integer: %d\n", intMultiItem->data.S32);
	}
	psFree(intIter);
    }
}
\end{verbatim}

The above program should produce, on being given no arguments:
\begin{verbatim}

Name of the program here

Usage: ./program INPUT OUTPUT

Optional arguments, with default values:
    -string    (SomeString)      Test string
    -bool      (FALSE)           Test bool
    -int       (1)               Test integer 1
               (2)               Test integer 2
               (3)               Test integer 3
    -float     (1.234567e+00)    Test float
\end{verbatim}

and called with: \code{./program -string whatevers -float 9.876543
-int 4 5 6 -bool input output}, then it should produce:

\begin{verbatim}
Success: input output
String: SomeString
Float: 1.234567
Boolean: 0
Integer: 1
Integer: 2
Integer: 3
\end{verbatim}


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

\pagebreak 
\section{Containers}

We require general data containers, so that associated values (e.g.\
the elements of an array) may be connected as a whole.  We require the
following types of containers:
\begin{itemize}
\item Arrays;
\item Doubly-linked lists;
\item Hashes;
\item Pixel lists;
\item Bit sets;
\item Metadata; and
\item Lookup tables.
\end{itemize}

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

\subsection{Arrays}

We require an order collection of unspecified data elements.  We
define \code{psArray} to carry such a collection:
%
\begin{datatype}
typedef struct {
    long n;                             ///< size of array 
    const long nalloc;                  ///< allocated data block
    psPtr *data;                        ///< pointer to data block
    void *lock;                         ///< Optional lock for thread safety
} psArray;
\end{datatype}
%
In this structure, the argument \code{n} is the length of the array
(the number of elements); \code{nalloc} is the number of elements
allocated ($nalloc \ge n$).  The allocated memory is pointed to by
\code{data}.  The structure is associated with a constructor and a
destructor:
%
\begin{prototype}
psArray *psArrayAlloc(long nalloc);
psArray *psArrayRealloc(psArray *array, long nalloc);
\end{prototype}
%
In these functions, \code{nalloc} is the number of elements to
allocate.  For \code{psArrayAlloc}, the value of \code{psArray.n} is
set to \code{nalloc}.  Users may choose to restrict the data range
after the \code{psArrayAlloc} function is called.  For
\code{psArrayRealloc}, if the value of \code{nalloc} is smaller than
the current value of \code{psArray.n}, then \code{psArray.n} is set to
\code{nalloc}, the array is adjusted down to match \code{nalloc}, and
the extra elements are dropped and freed if necesitated by the
reference counter.  If \code{nalloc} is larger than the current value
of \code{psArray.n}, \code{psArray.n} is left intact.  If the value of
\code{array} is \code{NULL}, then \code{psArrayRealloc} must return an
error.

\begin{prototype}
void psArrayElementsFree(psArray* array);
\end{prototype}
%
\code{psArrayElementsFree} shall free all elements on the
\code{array}.

\begin{prototype}
psArray *psArrayAdd(psArray *array, long delta, psPtr data);
\end{prototype}

This function adds a value to the end of an array.  If the current
length of the array (\code{psArray.n}) is at the limit of the
allocated space, additional space is allocated.  The value of
\code{delta} defines how many elements to add on each pass (if this
value is less than 1, 10 shall be used).

\begin{prototype}
bool psArrayRemove(psArray *array, const psPtr data);
\end{prototype}

This function removes all entries of \code{value} in the \code{array},
reducing the total number of elements of \code{array} as needed.
Returns \code{TRUE} if any elements were removed, otherwise
    const int x0, y0;                   ///< data region relative to parent 
\code{FALSE}.

\begin{prototype}
bool psArraySet(psArray *array, long position, psPtr data);
psPtr psArrayGet(const psArray *array, long position);
\end{prototype}

These accessor functions are provided as a convenience to the user.
\code{psArraySet} sets the value of the \code{in} array at the
specified \code{position} to \code{value}, returning \code{true} if
successful.  \code{psArrayGet} returns the value of the \code{in}
array at the specified \code{position}.  A negative \code{position}
means index from the end.

\begin{datatype}
typedef int (*psComparePtrFunc) (
    const void **a,                    ///< first comparison target
    const void **b                     ///< second comparison target
);
\end{datatype}

\begin{prototype}
psArray *psArraySort(psArray *array, psComparePtrFunc func);
\end{prototype}
An array may be sorted using \code{psArraySort}, which requires the
specification of a comparison function to specify how the objects on
the list should be sorted.  The motivation is primarily to be able to
iterate on a sorted list of keys from a hash.  The \code{array} is
sorted in-place.

\subsubsection{Comparison functions}

We specify several comparison functions as a convenience to the user.
The \code{psCompareTypePtr} functions are intended for sorting where
the container has a pointer to the value (e.g., \code{psList}).  The
functions return an integer less than, equal to, or greater than zero
if \code{a} is is less than, equal to, or greater than \code{b},
respectively (\code{qsort} man page).

\begin{prototype}
int psCompareS8Ptr(const void **a, const void **b);
int psCompareS16Ptr(const void **a, const void **b);
int psCompareS32Ptr(const void **a, const void **b);
int psCompareS64Ptr(const void **a, const void **b);
int psCompareU8Ptr(const void **a, const void **b);
int psCompareU16Ptr(const void **a, const void **b);
int psCompareU32Ptr(const void **a, const void **b);
int psCompareU64Ptr(const void **a, const void **b);
int psCompareF32Ptr(const void **a, const void **b);
int psCompareF64Ptr(const void **a, const void **b);
int psCompareDescendingS8Ptr(const void **a, const void **b);
int psCompareDescendingS16Ptr(const void **a, const void **b);
int psCompareDescendingS32Ptr(const void **a, const void **b);
int psCompareDescendingS64Ptr(const void **a, const void **b);
int psCompareDescendingU8Ptr(const void **a, const void **b);
int psCompareDescendingU16Ptr(const void **a, const void **b);
int psCompareDescendingU32Ptr(const void **a, const void **b);
int psCompareDescendingU64Ptr(const void **a, const void **b);
int psCompareDescendingF32Ptr(const void **a, const void **b);
int psCompareDescendingF64Ptr(const void **a, const void **b);
int psCompareS8(const void *a, const void *b);
int psCompareS16(const void *a, const void *b);
int psCompareS32(const void *a, const void *b);
int psCompareS64(const void *a, const void *b);
int psCompareU8(const void *a, const void *b);
int psCompareU16(const void *a, const void *b);
int psCompareU32(const void *a, const void *b);
int psCompareU64(const void *a, const void *b);
int psCompareF32(const void *a, const void *b);
int psCompareF64(const void *a, const void *b);
int psCompareDescendingS8(const void *a, const void *b);
int psCompareDescendingS16(const void *a, const void *b);
int psCompareDescendingS32(const void *a, const void *b);
int psCompareDescendingS64(const void *a, const void *b);
int psCompareDescendingU8(const void *a, const void *b);
int psCompareDescendingU16(const void *a, const void *b);
int psCompareDescendingU32(const void *a, const void *b);
int psCompareDescendingU64(const void *a, const void *b);
int psCompareDescendingF32(const void *a, const void *b);
int psCompareDescendingF64(const void *a, const void *b);
\end{prototype}

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

\subsection{Doubly-linked lists}
\label{sec:psList}

\PS{} shall support doubly linked lists through a type \code{psList}:
%
\begin{datatype}
typedef struct {
   long n;                              ///< number of elements on list
   psListElem *head;                    ///< first element on list (may be NULL)
   psListElem *tail;                    ///< last element on list (may be NULL)
   psArray *iterators;                  ///< array of psListIterator: iteration cursors
   void *lock;                          ///< Optional lock for thread safety
} psList;
\end{datatype}
%
The type \code{psList} represents the container of the list.  It has a
pointer to the first element in the linked list (\code{head}), a
pointer to the last element in the list (\code{tail}), an array of
iteration cursors, (\code{iterators}), and an entry to define the
number of elements in the list (\code{n}).

The elements of the list are defined by the type \code{psListElem}:
%
\begin{datatype}
typedef struct psListElem {
   struct psListElem *prev;            ///< previous link in list
   struct psListElem *next;            ///< next link in list
   psPtr data;                         ///< real data item
} psListElem;
\end{datatype}
%
which includes a pointer to the next element in the list
(\code{next}), the previous element in the list (\code{prev}), and a
\code{void} pointer to whatever data is represented by this list
element.    The following supporting functions are required:

\begin{prototype}
psList *psListAlloc(psPtr data);
\end{prototype}
Create a list.  This function may take a pointer to a data item, or it
may take \code{NULL}.  The allocator creates both the \code{psList}
and the first element in the list, pointed to by both
\code{psList.head} and \code{psList.tail}.  If the data entry is
\code{NULL}, then an empty list, with both pointers set to \code{NULL}
should be created.

The destructor function for \code{psList} must call \code{psFree} for
all the the data associated with the list.

All data items placed onto lists must have their reference counters
(section~\ref{secMemRefcounter}) incremented.  When elements are
removed from a list, they must have their reference counters
decremented.  The action of retrieving data from a list (with one of
the three \code{psListGet} functions) is considered ``borrowing'' the
reference, so no action is performed on the reference counter.

Iteration on the list shall be achieved by means of a list iterator
type:
\begin{datatype}
typedef struct {
    psList *list;                       ///< List iterator works on
    psListElem *cursor;                 ///< The current iterator cursor
    bool offEnd;                        ///< Is the iterator off the end?
    long index;                         ///< Index of iterator, to assist performance
    bool mutable;                       ///< Is it permissible to modify the list?
} psListIterator;
\end{datatype}
The \code{psListIterator} keeps track of which list element the
iterator \code{cursor} is currently pointing at.  \code{index} is the
index of the list iterator, which is used to assist performance when
using numerical locations.  The boolean member, \code{offEnd},
indicates whether the iterator has progressed off the end of the list
(i.e., beyond the last item).  The boolean \code{mutable} specifies
whether it is permissible to modify the list pointed to by the
iterator.  \code{psListAddBefore} and \code{psListAddAfter} are not
permitted to modify a list that is not \code{mutable} (i.e., only the
\code{psListGetAndIncrement} and \code{psListGetAndDecrement}
operations are permissible for a non-\code{mutable} list).

The corresponding constructor shall be:
\begin{prototype}
psListIterator *psListIteratorAlloc(psList *list, long location, bool mutable);
\end{prototype}
Here, \code{list} is the \code{psList} on which the iterator will
iterate, and \code{location} is the initial starting point, and may be
a numerical index or it may be one of the special values:
\code{PS_LIST_HEAD} or \code{PS_LIST_TAIL}, which are defined as 0 and
-1, respectively; a negative index is interpreted as relative to the
end of the list.  The boolean \code{mutable} specifies whether it is
permissible to modify the list pointed to by the iterator.

The destructor for \code{psListIterator} shall, after freeing the
\code{psListIterator}, also reorganise the \code{iter} array
(replacing the element being removed with the last element) and
resizing the array appropriately.

A list \code{iterator} shall be set to a specific \code{location} on
the list upon calling \code{psListIteratorSet}:
\begin{prototype}
bool psListIteratorSet(psListIterator *iterator, long location);
\end{prototype}
Again, the \code{location} may be a numerical index or it may be one
of the special values: \code{PS_LIST_HEAD} or \code{PS_LIST_TAIL},
which are defined as 0 and -1, respectively; a negative index is
interpreted as relative to the end of the list.  The function shall
return \code{true} if the reset was successful, or \code{false}
otherwise.

\begin{prototype}
bool psListAdd(psList *list, long location, psPtr data);
bool psListAddAfter(psListIterator *iterator, psPtr data);
bool psListAddBefore(psListIterator *iterator, psPtr data);
\end{prototype}
the first function, \code{psListAdd}, adds an entry to the \code{list}
and returns a boolean giving the success or failure of the
operation. The value of \code{location} may be a numerical index the
\code{data} is to inhabit (if \code{location} is greater than the
number of items on the list, then the function shall generate a
warning and add the \code{data} to the tail) or it may be one of the
special values: \code{PS_LIST_HEAD} or \code{PS_LIST_TAIL}, which are
defined as 0 and -1, respectively; a negative index is interpreted as
relative to the end of the list.  The other two functions,
\code{psListAddAfter} and \code{psListAddBefore} specify that the
\code{data} shall be added after or before (respectively) the current
cursor position of the \code{iterator}.

\begin{prototype}
psPtr psListGet(psList *list, long location);
psPtr psListGetAndIncrement(psListIterator *iterator);
psPtr psListGetAndDecrement(psListIterator *iterator);
\end{prototype}
A data item may be retrieved from the list with these functions.  The
first function, \code{psListGet} simply returns the value specified by
its \code{location}, which may be a numerical index or it may be one
of the special values: \code{PS_LIST_HEAD = 0} or \code{PS_LIST_TAIL =
-1}; negative indices are interpreted as relative to the end of the
list.  The other two functions, \code{psListGetAndIncrement} and
\code{psListGetAndDecrement} return the item under the iteration
cursor before advancing to the next or previous item, respectively.

In the event that the iteration cursor is at the tail of the list,
\code{psListGetAndIncrement} shall return the tail item and then set
the \code{cursor} to \code{NULL} and \code{offEnd} to \code{true}.  In
the event that the iteration cursor is at the head of the list,
\code{psListGetAndDecrement} shall return the head item and then set
the \code{cursor} to \code{NULL} (and \code{offEnd} should already be
\code{false}).  In the event that the iteration \code{cursor} is
\code{NULL}, \code{psListGetAndIncrement} and
\code{psListGetAndDecrement} shall return \code{NULL}, and advance the
iteration \code{cursor} only if the intended direction places the
cursor back on the list; otherwise a warning shall be generated, and
no change shall be made.  If \code{psListGetAndDecrement} was called
with \code{offEnd} as \code{true}, then \code{offEnd} shall also be
toggled back to \code{false} to indicate that the \code{cursor} is no
longer off the end of the list.

\begin{prototype}
bool psListRemove(psList *list, long location)
bool psListRemoveData(psList *list, psPtr data);
\end{prototype}
A data item may be removed from the list with these functions.  For
\code{psListRemove}, the value of \code{location} may be the numerical
index or it may be one of the special values: \code{PS_LIST_HEAD} or
\code{PS_LIST_TAIL}, which are defined as 0 and -1, respectively; a
negative index is interpreted as relative to the end of the list.  For
\code{psListRemoveData}, the data item to be removed is identified by
matching the pointer value with \code{psPtr data}.  The functions
return a value of \code{true} if the operation was successful, and
\code{false} otherwise.  In both cases, if any iterators are currently
pointing at the item to be removed, the item shall be removed and
those iterators pointing at it shall be moved to the next, and the
function shall return \code{true}.  If the item to be removed is not
on the list, an error shall be generated and the function shall return
\code{false}.

\begin{prototype}
psArray *psListToArray(const psList *list);
psList  *psArrayToList(const psArray *array);
\end{prototype}
These two functions are available to convert between the
\code{psList} and \code{psArray} containers.  These functions do not
free the elements or destroy the input collection.  Rather, they
increment the reference counter for each of the elements.

\begin{prototype}
psList *psListSort(psList *list, psComparePtrFunc func);
\end{prototype}
A list may be sorted using \code{psListSort}, which requires the
specification of a comparison function to specify how the objects on
the list should be sorted.  The motivation is primarily to be able to
iterate on a sorted list of keys from a hash.  The \code{list} is
sorted in-place.

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

\subsection{Hash Tables}
\hlabel{psHash}

Hash tables are critical for quick retrieval of text-based data.  The
concept is as follows: Given a large collection of text strings, it is
inefficient to search for a particular entry by performing a basic
string comparison on all entries until a match is found.  Even if a
single list is sorted, we will still spend a substantial amount of
time iterating across the entries in the list.  In a hash table, we
define an operation, the hash function, which uses the bytes of the
string to construct a numerical value, the hash value.  The hash value
is defined to have a limited range of $N$ values.  The hash table
consists of $N$ buckets, each of which contains a list of the strings
whose hash value corresponds to the bucket number.  Searching for a
specific string involves calculating the hash value for the string,
going to the appropriate bucket, and searching through the
corresponding list until the string is matched.  

For PSLib, we define a hash table and hash buckets as follows:
\footnote{ We choose not to use the POSIX function \code{hcreate},
\code{hdestroy}, and \code{hsearch} as they only support a single hash
table at any one time.}
%
\begin{datatype}
typedef struct {
    long n;                             ///< number of buckets
    psHashBucket **buckets;             ///< the buckets themselves
    void *lock;                         ///< Optional lock for thread safety
} psHash;
\end{datatype}
%
where \code{n} is the number of buckets defined for the hash functions, and
\code{buckets} is an array of pointers to the individual buckets, each of which
is defined by:
%
\begin{datatype}
typedef struct psHashBucket {
    char *key;                          ///< key for this item of data
    psPtr data;                         ///< the data itself
    struct psHashBucket *next;          ///< list of other possible keys
} psHashBucket;
\end{datatype}
where each bucket contains the value of the \code{key}, a pointer to
the \code{data}, and a pointer to the \code{next} list entry in the
bucket (in the event that two or more keys have the same hash value).

A hash table is created with the following function:
\begin{prototype}
psHash *psHashAlloc(long nalloc);
\end{prototype}
which allocates the space for the hash table, creating and
initializing \code{n} hash buckets.

The destructor for \code{psHash} must free all data associated with a complete hash table.

A data item may be added to the hash table with the function:
\begin{prototype}
bool psHashAdd(psHash *hash, const char *key, psPtr data);
\end{prototype}
In this function, the value of the string \code{key} is used to
construct the hash value, find the appropriate bucket set, and add the
new element \code{data} to the list.  An existing element with the same
value of \code{key} is destroyed using its registered destructor
(\code{psMemBlock}). The return value of the function is a boolean
defining the success or failure of the operation.

The data associated with a given key may be found with the function:
\begin{prototype}
psPtr psHashLookup(const psHash *hash, const char *key);
\end{prototype}
which returns the data value pointed to by the element associated with
\code{key}, or the value \code{NULL} if no match is found.  Similarly,
a specific key may be removed (deleted) with the function:
\begin{prototype}
bool psHashRemove(psHash *hash, const char *key);
\end{prototype}
The function returns a value of \code{true} if the operation was
successfull, and \code{false} otherwise.

The function
\begin{prototype}
psList *psHashKeyList(const psHash *hash);
\end{prototype}
returns the complete list of defined keys associated with the
\code{psHash} table as a linked list.

\begin{prototype}
psArray *psHashToArray(const psHash *hash);
\end{prototype}
This function places the data in a \code{psHash} into a \code{psArray}
container.  This function does not free the elements or destroy the
input collection.  Rather, it increments the reference counter for
each of the elements.  The resulting array does not have any
information about the has key values, and the order is not
significant.

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

\subsection{Pixel Lists}

Usually an image mask is the best way to carry information about what
pixels mean what.  However, in the case where the number of pixels in
which we are interested is limited, it is more efficient to simply
carry a list of pixels.  An example of this is in the image
combination code, where we want to perform an operation on a
relatively small fraction of pixels, and it is inefficient to go
through an entire mask image checking each pixel.

\begin{datatype}
typedef struct {
    float x;                    // x coordinate
    float y;                    // y coordinate
} psPixelCoord;

typedef struct {
    long n;                     // Number in use
    const long nalloc;          // Number allocated
    psPixelCoord *data;         // The pixel coordinates
    void *lock;                 // Lock for thread safety
} psPixels;
\end{datatype}

\begin{prototype}
psPixels *psPixelsAlloc(long nalloc);
psPixels *psPixelsRealloc(psPixels *pixels, long nalloc);
\end{prototype}

\code{psPixelsAlloc} and \code{psPixelsRealloc} provide dynamic
allocation and reallocation in a manner analogous to those provided
by \code{psVectorAlloc} and \code{psVectorRealloc}.

\begin{prototype}
psPixels *psPixelsCopy(psPixels *out, const psPixels *pixels);
psPixels *psPixelsConcatenate(psPixels *out, const psPixels *pixels);
\end{prototype}

\code{psPixelsCopy} shall copy the contents of \code{pixels} to the
\code{out}.  In the event that \code{out} is \code{NULL}, a new
\code{psPixels} shall be allocated, and the contents of \code{pixels}
simply copied in.  If \code{pixels} is \code{NULL}, the function shall
generate an error and return \code{NULL}.

\code{psPixelsConcatenate} shall concatenate the \code{pixels} onto
\code{out}.  In the event that \code{out} is \code{NULL}, the function
performs a \code{psPixelsCopy}, returning the copy.  If \code{pixels}
is \code{NULL}, the function shall generate an error and return
\code{NULL}.  The function shall take care to ensure that there are no
duplicate pixels in \code{out} (since the order in which the pixels
are stored is not important, the values may be sorted, allowing the
use of a faster algorithm than a linear scan).

\begin{prototype}
bool psPixelsSet(psPixels *pixels, long position, psPixelCoord value);
psPixelCoord psPixelsGet(const psPixels *pixels, long position);
\end{prototype}

These accessor functions are provided as a convenience to the user.
\code{psPixelsSet} sets the value of the \code{pixels} array at the
specified \code{position} to \code{value} (a \code{psPixelCoord}
passed by value), returning \code{true} if successful.
\code{psPixelsGet} returns the value of the \code{pixels} array at the
specified \code{position}.  A negative \code{position} means index
from the end.

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

\subsection{BitSets}

BitSets are required in order to turn options on and off.  We require
the capability to have a bitset of arbitrary length (i.e., not limited
by the length of a \code{long}, say).  The \code{psBitSet} structure
is defined below.  Note that the entry \code{bits} is an array of type
\code{char} storing the bits as bits of each byte in the array, with 8
bits available for each byte in the array.  Also note that the
constructor is passed the number of required bits, which implies that
\code{ceil(n/8)} bytes must be allocated.  The bitset structure is
define by:
\begin{datatype}
typedef struct {
    long n;                             ///< Number of chars that form the bitset
    psU8 *bits;                         ///< The bits
    void *lock;                         ///< Optional lock for thread safety
} psBitSet;
\end{datatype}

We also require the corresponding constructor and destructor:
\begin{prototype}
psBitSet *psBitSetAlloc(long nalloc);
\end{prototype}
where \code{n} is the requested number of bits.

Several basic operations on bitsets are required:
\begin{itemize}
\item Set a bit;
\item Check if a bit is set; and
\item \code{OR}, \code{AND} and \code{XOR} two bitsets.
\item \code{NOT} a bitset.
\end{itemize}
The corresponding APIs are defined below:

\begin{prototype}
psBitSet *psBitSetSet(psBitSet *bitSet, long bit);
psBitSet* psBitSetClear(psBitSet *bitSet, long bit);
psBitSet *psBitSetOp(psBitSet *outBitSet, const psBitSet *inBitSet1, const char *operator, const psBitSet *inBitSet2);
psBitSet *psBitSetNot(psBitSet *outBitSet, const psBitSet *inBitSet);
bool psBitSetTest(const psBitSet *bitSet, long bit);
psString psBitSetToString(const psBitSet* bitSet);
\end{prototype}

\code{psBitSetSet} sets the specified \code{bit} in the
\code{psBitSet}, and returns the updated bitset.  The input bitset
will be modified.

\code{psBitSetClear} clears the specified \code{bit} in the \code{bitSet}
and returns the updated bitset.  The input bitset will be modified.

\code{psBitSetOp} returns the \code{psBitSet} that is the result of
performing the specified \code{operator} (one of \code{"AND"},
\code{"OR"}, or \code{"XOR"}) on \code{inBitSet1} and \code{inBitSet2}.
If the output bitset \code{outBitSet} is \code{NULL}, it is created by
the function.

\code{psBitSetNot} applies a unary \code{NOT} to a bitset, placing the
answer in the bitset \code{out}, or creating a new bitset if
\code{out} is \code{NULL}.

\code{psBitSetTest} returns a true value if the specified \code{bit}
is set; otherwise, it returns a false value.

Finally, \code{psBitSetToString} returns a string representation of
the specified \code{bits}.

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

\subsection{Metadata}
\label{sec:metadata}

\subsubsection{Conceptual Overview}

Within PSLib, we provide a data structure to carry metadata and
mechanisms to manipulate the metadata.  Metadata is a general concept
that requires some discussion.  In any data analysis task, the
ensemble of all possible data may be divided into two or three
classes: there is the specific data of interest, there is data which
is related or critical but not the primary data of interest, and there
is all of the other data which may or may not be interesting.  For
example, consider a simple 2D image obtained of a galaxy from a CCD
camera on a telescope.  If you want to study the galaxy, the specific
data of interest is the collection of pixels.  There are a variety of
other pieces of data which are closely related and crucial to
understanding the data in those pixels, such as the dimensions of the
image, the coordinate system, the time of the image, the exposure
time, and so forth.  Other data may be known which may be less
critical to understanding the image, but which may be interesting or
desired at a later date.  For example, the observer who took the
image, the filter manufacturer, the humidity at the telescope, etc.

Formally, all of the related data which describe the principal data of
interest are metadata.  Note that which piece is the metadata and
which is the data may depend on the context.  If you are examining the
pixels in an image, the coordinate and flux of an object may be part
of the metadata.  However, if you are analyzing a collection of
objects extracted from an image, you may consider then pixel data
simply part of the metadata associated with the list of objects.  

There are various ways to handle metadata vs data within a programming
environment.  In C, it is convenient to use structures to group
associated data together.  One possibility is to define the metadata
as part of the associated data structure.  For example, the image data
structure would have elements for all possible associated measurement.
This approach is both cumbersome (because of the large number metadata
types), impractical (because the full range of necessary metadata is
difficult to know in advance) and inflexible (because any change in
the collection of metadata requires addition of new structure elements
and recompilation).  

An alternative is to place the metadata in a generic container and use
lookup mechanisms to extract the appropriate metadata when needed.  An
example of this is would be a text-based FITS header for an image read
into a flat text buffer.  In this implementation, metadata lookup
functions could return the current value of, for example, NAXIS1 (the
number of columns of the image) by scanning through the header buffer.
This method has the benefits of flexibility and simplicity of
programming interface, but it has the disadvantage that all metadata
is accessed though this lookup mechanism.  This may make the code less
readable and it may slow down the access.  

PSLib implements an intermediate solution to this problem.  We specify
a flexible, generic metadata container and access methods.  Data types
which require association with a general collection of metadata should
include an entry of this metadata type.  However, a subset of metadata
concepts which are basic and frequently required may be placed in the
coded structure elements.  This approach allows the code to refer to
the basic metadata concepts as part of the data structure (ie,
\code{image.nx}), but also allows us to provide access to any
arbitrary metadata which may be generated.  As a practical matter, the
choice of which entries are only in the metadata and which are part of
the explicit structure elements is rather subjective.  Any data
elements which are frequently used should be put in the structure;
those which are only infrequently needed should be left in the generic
metadata.

There are some points of caution which must be noted.  Any
manipulation of the data should be reflected in the metadata where
appropriate.  This is always an issue of concern.  For example,
consider an image of dimensions \code{nx, ny}.  If a function extracts
a subraster, it must change the values of \code{nx, ny} to match the
new dimensions.  What should it do to the corresponding metadata?
Clearly, it should change the corresponding value which defines
\code{nX, nY}.  However, it is not quite so simple: there may be other
metadata values which depend on those values.  These must also be
changed appropriately.  What if the metadata element points to a
copy of the metadata which may be shared by other representations of
the image?  These must be treated differently because the change would
invalidate those other references.  Care must be taken, therefore,
when writing functions which operate on the data to consider all of
the relevant metadata entries which must also be updated. 

A related issue is the definition of metadata names.  Entries in a
structure have the advantage of being hardwired: every instance of
that structure will have the same name for the same entry.  This is
not necessarily the case with a more flexible metadata container.  The
image exposure time is a notorious example in astronomy.  Different
observatories use different header keywords (ie, metadata names) for
the same concept of the exposure time (\code{EXPTIME},
\code{EXPOSURE}, \code{OPENTIME}, \code{INTTIME}, etc).  Any system
which operates on these metadata needs to address the issue of
identifying these names.  This issue seems like an argument for
hardwiring metadata in the structure, but in fact it does not present
such a strong case.  If the metadata are hardwired, some function will
still have to know how to interpret the various names to populate the
structure.  The concept can still be localized with generic metadata
containers by including abstract metadata names within the code which
are tied to the various implementations-specific metadata names.

\subsubsection{Metadata Representation}

\begin{figure}
\psfig{file=Metadata,width=6.5in}
\caption{Metadata Structures\label{fig:metadata}}
\end{figure}

This section addresses the question of how \PS{} metadata should be
represented in memory, not how it should be represented on disk.

We define an item of metadata with the following structure:
\filbreak
\begin{datatype}
typedef struct {
    const psS32 id;                     ///< unique ID for this item
    psString name;                      ///< Name of item
    psDataType type;                    ///< type of this item
    union {
        psBool B;                       ///< boolean value
        psS8  S8;                       ///< integer data
        psS16 S16;                      ///< integer data
        psS32 S32;                      ///< integer data
        psU8  U8;                       ///< integer data
        psU16 U16;                      ///< integer data
        psU32 U32;                      ///< integer data
        psF32 F32;                      ///< floating-point data
        psF64 F64;                      ///< double-precision data
        psList *list;                   ///< psList entry
        psMetadata *md;                 ///< psMetadata entry
        psPtr V;                        ///< other type
    } data;                             ///< value of metadata
    psString comment;                   ///< optional comment ("", not NULL)
} psMetadataItem;
\end{datatype}

The \code{id} is a unique identifier for this item of metadata;
experience shows that such tags are useful.  The entry \code{name}
specifies the name of the metadata item.  The value of the metadata is
given by the union \code{data}, and may be of type \code{psU8},
\code{psU16}, \code{psU32}, \code{psS8}, \code{psS16}, \code{psS32},
\code{psF32}, \code{psF64}, or an arbitrary rich structure pointed at
by the \code{void} pointer \code{V}.  A character string comment
associated with this metadata item may be stored in the element
\code{comment}. The \code{type} entry specifies how to interpret the
type of the data being represented, given by the enumerated type
\code{psDataType}.

Note that the \code{name} and \code{comment} must be allocated by the
constructor using, e.g., \code{psStringCopy}.

A collection of metadata is represented by the \code{psMetadata} structure:
\begin{datatype}
typedef struct {
    psList *list;                       ///< list of psMetadataItem
    psHash *hash;                       ///< hash table of the same metadata
    void *lock;                         ///< Optional lock for thread safety
} psMetadata;
\end{datatype}
The type \code{psMetadata} is a container class for metadata. Note
that there are in fact \emph{two} representations of the metadata
(each \code{psMetadataItem} appears on both).  The first
representation employs a doubly-linked list that allows the order of
the metadata to be preserved (e.g., if FITS headers are read in a
particular order, they should be written in the same order).  The
second representation employs a hash table which allows fast look-up
given a specific metadata keyword.

Certain metadata names (such as the FITS keywords \code{COMMENT} and
\code{HISTORY} in a FITS header) may be repeated with different
values.  In such a case, the \code{psMetadata.list} structure contains
the entries in their original sequence with duplicate keys.  The
\code{psMetadata.hash} entries, which are required to have unique
keys, would have a single entry with the keyword of the repeated key,
with the value of \code{psDataType} set to \code{PS_DATA_METADATA_MULTI},
and the \code{psMetadataItem.data} element pointing to a \code{psList}
containing the actual entries.  If \code{psMetadataItemAlloc} is
called with the type set to \code{PS_DATA_METADATA_MULTI}, such a repeated key
is created.  In this case, the data value passed to
\code{psMetadataItemAlloc} (the quantity in ellipsis) must be
\code{NULL}.  An empty \code{psMetadataItem} with the given keyword is
created to hold future entries of that keyword.

As a convenience to the user, the following type-specific functions are
also defined:
\begin{prototype}
psMetadataItem* psMetadataItemAllocStr(const char* name, const char* comment, const char* value);
psMetadataItem* psMetadataItemAllocF32(const char* name, const char* comment, psF32 value);
psMetadataItem* psMetadataItemAllocF64(const char* name, const char* comment, psF64 value);
psMetadataItem* psMetadataItemAllocS32(const char* name, const char* comment, psS32 value);
psMetadataItem* psMetadataItemAllocBool(const char* name, const char* comment, bool value);
psMetadataItem* psMetadataItemAllocPtr(const char* name, psDataType type, const char* comment, psPtr value);
\end{prototype}

\subsubsection{Metadata APIs}

\begin{prototype}
psMetadata *psMetadataAlloc(void);
\end{prototype}

The constructor for the collection of metadata, \code{psMetadata},
simply returns an empty metadata container (employing the constructors
for the doubly-linked list and hash table).  The destructor needs to
free each of the \code{psMetadataItem}s.

\begin{prototype}
psMetadataItem *psMetadataItemAlloc(const char *name, psDataType type, const char *comment, ...);
psMetadataItem *psMetadataItemAllocV(const char *name, psDataType type, const char *comment, va_list list);
\end{prototype}

The allocator for \code{psMetadataItem} returns a full
\code{psMetadataItem} ready for insertion into the \code{psMetadata}.
The \code{name} entry specifies the name to use for this metadata
item, and may include \code{sprintf}-type formating codes.  The
\code{comment} entry is a fixed string which is used for the comment
associated with this metadata item.  The metadata data and the
arguments to the \code{name} formatting codes are passed, in that
order (metadata pointer first), to \code{psMetadataItemAlloc} as
arguments following the comment string.  The data must be a pointer
for any data types which are stored in the element \code{data.void},
while other data types are passed as numeric values.  The argument
list must be interpreted appropriately by the \code{va_list} operators
in the function.

\begin{prototype}
bool psMetadataAddItem(psMetadata *md, const psMetadataItem *item, int location, psS32 flags);
bool psMetadataAdd(psMetadata *md, long location, const char *name, int format, const char *comment, ...);
bool psMetadataAddV(psMetadata *md, long location, const char *name, int format, const char *comment,
                    va_list list);
\end{prototype}

Items may be added to the metadata in one of two ways --- firstly, an
item may be added by appending a \code{psMetadataItem} which has
already been created; and secondly by directly providing the data to
be appended.  In both cases, the return value defines the success
(\code{true}) or failure of the operation.  The second function,
\code{psMetadataAdd} takes a pointer or value which is interpreted by
the function using variadic argument interpretation.  The third
version is the \code{va_list} version of the second function.  All
three functions take a parameter, \code{location}, which specifies
where in the list to place the element, following the conventions for
the \code{psList}.  The entry \code{mode} for \code{psMetadataAddItem}
is a bit mask constructed by OR-ing the allowed option flags (eg,
\code{PS_DATA_REPLACE}) which specify minor variations on the
behavior.  The \code{format} entry, which specifies both the metadata
type and the optional flags, is constructed by bit-wise OR-ing the
appropriate \code{psDataType} and allowed option flags.  Care
should be taken not to leak memory when appending an item for which
the key already exists in the metadata (and is not
\code{PS_DATA_METADATA_MULTI}).
%

\begin{datatype}
typedef enum {                          ///< option flags for psMetadata functions
    PS_META_DEFAULT         = 0,        ///< default behavior (0x0000) for use in mode above
    PS_META_REPLACE         = 0x1000000, ///< allow entry to be replaced
    PS_META_DUPLICATE_OK    = 0x2000000, ///< allow duplicate entries
    PS_META_NULL            = 0x4000000 ///< psMetadataItem.data is a NULL value
} psMetadataFlags;
\end{datatype}

The functions above take option flags which modify the behavior when
metadata items are added to the metadata list.  These flags must be
bit-exclusive of those used above for the \code{psDataTypes}.  The
flags have the following meanings: 

\code{PS_META_DEFAULT}: This is the zero bit mask, to allow the
default behavior for \code{psMetadataAddItem} above.  If this is OR-ed
with a \code{psDataType}, the result is as if no OR-ing took
place.

\code{PS_META_REPLACE}: Replace an existing, unique entry. If the
given metadata item exists in the metadata collection, and is not of
type \code{PS_META_MULTI}, then the item replaces the existing entry.

\code{PS_META_DUPLICATE_OK}: Allow the new metadata item key to be a
duplicate (ie, \code{PS_DATA_METADATA_MULTI}).  If an existing item
with the same key is already \code{PS_DATA_METADATA_MULTI}, the new
item is added to the \code{PS_DATA_METADATA_MULTI} list.  If the
existing item is not \code{PS_DATA_METADATA_MULTI}, a
\code{PS_DATA_METADATA_MULTI} list is created to contain both the
existing item and the new item.  The original entry's location on the
psMetadata.list must be maintained.

\code{PS_META_NULL}: Indicates that \code{psMetadataItem.data} should
be ignored and that the the current value is ``NULL'' or undefined.
The \code{psMetadataItem} must have a proper \code{type} set and it's
\code{data} field shall have a valid value.  e.g. A \code{type} of
\code{PS_DATA_STR} would require that its \code{data} is set to
\code{NULL}.

There are several of cases to handle for duplication of an existing
key by a new key, some identified above.  The following situations
must also be handled:

If the new key already exists, but is not
\code{PS_DATA_METADATA_MULTI}, and the new item is not flagged as
either \code{PS_META_DUPLICATE_OK} or \code{PS_META_REPLACE}, an error
is raised.

If the new key already exists, and the existing item is
\code{PS_DATA_METADATA_MULTI}, the new item is added to the MULTI
list.  Note that if the new item is also of type
\code{PS_DATA_METADATA_MULTI}, no action is taken, but a successful
exit status is returned (the action of adding a
\code{PS_DATA_METADATA_MULTI} item to the metadata is equivalent to
setting that key to be tagged as \code{PS_DATA_METADATA_MULTI}.  If it
is {\em already} \code{PS_DATA_METADATA_MULTI}, this effect has
already been achieved).

An example of code to use these metadata APIs to generate the
structure seen in Figure~\ref{fig:metadata} is given below.

\begin{verbatim}
md = psMetadataAlloc();

psMetadataAdd(md, PS_LIST_TAIL, "SIMPLE",   PS_DATA_BOOL, "basic fits",            TRUE);
psMetadataAdd(md, PS_LIST_TAIL, "BLANK",    PS_DATA_S32,  "invalid pixel data",    -32768);
psMetadataAdd(md, PS_LIST_TAIL, "DATE-OBS", PS_DATA_STR,  "observing date UT", "   2004-6-16");
psMetadataAdd(md, PS_LIST_TAIL, "COMMENT",  PS_DATA_LIST, "head of comment block", NULL);
psMetadataAdd(md, PS_LIST_TAIL, "COMMENT",  PS_DATA_STR,  "",                      "DATA");
psMetadataAdd(md, PS_LIST_TAIL, "COMMENT",  PS_DATA_STR,  "",                      "PARAMS"); 
psMetadataAdd(md, PS_LIST_TAIL, "EXPTIME",  PS_DATA_F32,  "exposure time (sec)",   1.05);
psMetadataAdd(md, PS_LIST_TAIL, "COMMENT",  PS_DATA_STR,  "",                      "FOO");

cell = psMetadataAlloc();
psMetadataAdd(cell, PS_LIST_TAIL, "EXTNAME",  PS_DATA_STR,  "",                    "CCD00");
psMetadataAdd(cell, PS_LIST_TAIL, "BIASNAME", PS_DATA_STR,  "",                    "BSEC-00");
psMetadataAdd(cell, PS_LIST_TAIL, "CHIP",     PS_DATA_STR,  "",                    "CHIP.00");
psMetadataAdd(md,   PS_LIST_TAIL, "CELL.00",  PS_DATA_META, "",                    cell);

cell = psMetadataAlloc();
psMetadataAdd(cell, PS_LIST_TAIL, "EXTNAME",  PS_DATA_STR,  "",                    "CCD01");
psMetadataAdd(cell, PS_LIST_TAIL, "BIASNAME", PS_DATA_STR,  "",                    "BSEC-01");
psMetadataAdd(cell, PS_LIST_TAIL, "CHIP",     PS_DATA_STR,  "",                    "CHIP.01");
psMetadataAdd(md,   PS_LIST_TAIL, "CELL.01",  PS_DATA_META, "",                    cell);
\end{verbatim}

The following code shows how to use the APIs to replace one of these values:
\begin{verbatim}
psMetadataAdd(md, PS_LIST_TAIL, "EXPTIME",  PS_DATA_F32 | PS_REPLACE,  "new exposure time (sec)",   2.05);
\end{verbatim}

As a convenience to the user, the following type-specific functions
are specified:
\begin{prototype}
bool psMetadataAddStr(psMetadata* md, long location, const char* name, int format,
                      const char* comment, const char* value);
bool psMetadataAddS32(psMetadata* md, long location, const char* name, int format,
                      const char* comment, psS32 value);
bool psMetadataAddF32(psMetadata* md, long location, const char* name, int format,
                      const char* comment, psF32 value);
bool psMetadataAddF64(psMetadata* md, long location, const char* name, int format,
                      const char* comment, psF64 value);
bool psMetadataAddBool(psMetadata* md, long location, const char* name, int format,
                       const char* comment, bool value);
bool psMetadataAddPtr(psMetadata* md, long location, const char* name, psDataType type,
                      const char* comment, psPtr value);
\end{prototype}


Items may be removed from the metadata by specifying a key or a
location in the list.  For \code{psMetadataRemoveKey}, if the key
matches a metadata item, the item is removed from the metadata and
\code{true} is returned; otherwise, \code{false} is returned.  If the
key is not unique, then \emph{all} items corresponding to the key are
removed, and \code{true} is returned.  For
\code{psMetadataRemoveIndex}, the metadata item at the specified
\code{location} is removed, if valid, and \code{true} is returned;
otherwise the function returns \code{false}.
%
\begin{prototype}
bool psMetadataRemoveKey(psMetadata *md, const char *key);
bool psMetadataRemoveIndex(psMetadata *md, int location);
\end{prototype}

Items may be found within the metadata by providing a key.  In the
event that the key is non-unique, the first item is returned.
\begin{prototype}
psMetadataItem *psMetadataLookup(const psMetadata *md, const char *key);
\end{prototype}

Several utility functions are provided for simple cases.  These
functions perform the effort of casting the data to the appropriate
type.  The numerical functions shall return 0.0 if their key is not
found.  If the pointer value of \code{status} is not \code{NULL}, it
is set to reflect the success or failure of the lookup.
\begin{prototype}
psPtr psMetadataLookupStr(bool *status, const psMetadata *md, const char *key);
psS32 psMetadataLookupS32(bool *status, const psMetadata *md, const char *key);
psF32 psMetadataLookupF32(bool *status, const psMetadata *md, const char *key);
psF64 psMetadataLookupF64(bool *status, const psMetadata *md, const char *key);
bool psMetadataLookupBool(bool *status, const psMetadata *md, const char *key);
psPtr psMetadataLookupPtr(bool *status, const psMetadata *md, const char *key);
\end{prototype}

Items may be retrieved from the metadata by their entry position.  The
value of which specifies the desired entry in the fashion of
\code{psList}.
\begin{prototype}
psMetadataItem *psMetadataGet(const psMetadata *md, int location);
\end{prototype}

The metadata list component may be iterated over by using a
\code{psMetadataIterator} in a fashion equivalent to the
\code{psListIterator}:
\begin{datatype}
typedef struct {
    psListIterator* iter;              ///< iterator for the psMetadata's psList
    regex_t* regex;                     ///< the subsetting regular expression
} psMetadataIterator;
\end{datatype}

The iterator may be set to a location in the \code{psMetadata} list,
and the user may get the previous or next item in the list relative to
that location.  The iterators may be used to return the next key
matching a POSIX \code{regex}, e.g., if the user only wants to iterate
through \code{IPP.machines.sky} and doesn't want to bother with
\code{IPP.machines.detector}.  The iterator should iterate over every
item in the metadata list, even those that are contained in a
\code{PS_DATA_LIST}.  The value \code{iterator} specifies the iterator
to be used.  In setting the iterator, the position of the iterator is
defined by \code{location}, which follows the conventions of the
\code{psList} iterators.
\begin{prototype}
psMetadataIterator *psMetadataIteratorAlloc(psMetadata *md, long location, const char *regex);
bool psMetadataIteratorSet(psMetadataIterator *iterator, long location);
psMetadataItem *psMetadataGetAndIncrement(psMetadataIterator *iterator);
psMetadataItem *psMetadataGetAndDecrement(psMetadataIterator *iterator);
\end{prototype}

Metadata items may be printed to an open file descriptor based on a
provided format.  The format string is an sprintf format statement
with exactly one \% formatting command.  If the metadata item type is
a numeric type, this formatting command must also be numeric, and type
conversion performed to the value to match the format type.  If the
metadata item type is a string, the formatting command must also be
for a string (\%s type of command).  If the metadata type is any other
data type, printing is not allowed.
\begin{prototype}
bool psMetadataItemPrint(FILE *fd, const char *format, const psMetadataItem *item);
\end{prototype}

\subsubsection{Configuration files}
\label{sec:configspec}

It will be necessary for the \PS{} system, in order to load
pre-defined settings, to parse a configuration file into a
\code{psMetadata} structure.  This shall be performed by the
function \code{psMetadataConfigParse}, as described below.

\begin{prototype}
psMetadata *psMetadataConfigParse(psMetadata *md, unsigned int *nFail, const char *filename, bool overwrite);
\end{prototype}

Given a metadata container, \code{md}, and the name of a configuration
file, \code{filename}, \code{psMetadataConfigParse} shall parse the
configuration file, placing the contained key/type/value/comment quads
into the metadata, and returning a pointer to the metadata structure.
The number of lines that failed to parse is returned in \code{nFail}.
Multiple specifications of a key that haven't been declared (see
below) are overwritten if and only if \code{overwrite} is \code{true}.
If the metadata container is \code{NULL}, it shall be allocated.  

On error, the function shall return \code{NULL}.

It is also useful to be able to convert a \code{psMetadata} structure into the
Configuration File format for debugging purposes and to enable persistent
configuration.

\begin{prototype}
psString psMetadataConfigFormat(psMetadata *md);
bool psMetadataConfigWrite(psMetadata *md, const char *filename);
\end{prototype}

The \code{psMetadataConfigFormat} function converts a \code{psMetadata}
structure (including any nested \code{psMetadata}) into a Configuration File
formatted string.  A \code{NULL} shall be returned on error.  The
\code{psMetadataConfigWrite} behaves the same as \code{psMetadataConfigFormat}
except that the string is written out to \code{filename}.  \code{false} is
returned on failure.

\paragraph{Comments}

The configuration file shall consist of plain text with
key/type/value/comment quads on separate lines.  Blank lines,
including those consisting solely of whitespace (both spaces and
tabs), shall be ignored, as shall lines that commence with the comment
character (a hash mark, \code{#}), either immediately at the start of
the line, or preceded by whitespace.  The key/type/value/comment quads
shall all lie on a single line, separated by whitespace.

The key shall be first, possibly preceded on the line by whitespace
which should not form part of the key.

\paragraph{NULL values}

The ``value'' of a quad may be declare to be undefined with the \code{NULL}
keyword.  \code{NULL} is allowed to co-exist with a ``comment'' and may be
surrounded by whitespace.  Any non-whitespace character will cause of the
``value'' to be interpreted as a string.

\begin{verbatim}
foo     STR     NULL    # string with a NULL value
bar     STR     NULL a  # string with a value of "NULL a"
\end{verbatim}

\paragraph{Types}
\subparagraph{Scalar \& Vector}

Next, to assist the casting of the value, shall be a string identifying the
type of the value, which shall correspond to one of the simple types supported
in \code{psMetadata}: \code{STRING,BOOL,S32,F32,F64}; \code{STR} may be used to
abbreviate \code{STRING}; valid time types are \code{UTC,UT1,TAI,TT}.

\tbd{May, in the future, require more types, including U8,S16,C64,
which will also necessitate updating the definition of psMetadata.}

The value shall follow the type: strings may consist of multiple words, and
shall have all leading and trailing whitespace removed; booleans shall simply
be either \code{T} or \code{F}.  Time type values will be in the ISO8601
compatible format of "YYYY-MM-DDTHH:MM:SS,sZ".  When parsed, time types shall
be represented as a \code{psTime} object.

Following the value may be an optional comment, preceded by a comment
character (a hash mark, \code{#}), which in the case of a string
value, serves to mark the end of the value, and for other types serves
to identify the comment to the reader.  Only one comment character may
be present on any single line (i.e., neither strings nor comments are
permitted to contain the comment character).  The comment may consist
of multiple words, and shall have leading and trailing whitespace
removed.

One wrinkle is the specification of vectors.  Keys for which the value
is to be parsed as a vector shall be preceded immediately by a
``vector symbol'', which we choose to be the ``at'' sign, \code{@}.
In this case, the type shall be interpreted as the type for the
vector, which may be any of the signed or unsigned integer or floating
point types (\code{U8,U16,U32,U64,S8,S16,S32,S32,S64,F32,F64}) but not
the complex floating point types; and the value shall consist of
multiple numbers, separated either by a comma or whitespace.  These
values shall populate a \code{psVector} of the appropriate type in the
order in which they appear in the configuration file.

\tbd{May add complex types, likely to be specified with values such as
  1.23+4.56i in the future.}

\tbd{May add null, Not-a-Number (NaN), de-normalized, underflow, overflow,
and/or +/-infinity values for selected types.}

\subparagraph{MULTI}

An additional hurdle is the specification of keys that may be non-unique (such
as the \code{COMMENT} keyword in a FITS header).  These keys shall be specified
in the configuration file as non-unique with a \code{MULTI} declaration.  In
the form \code{[keyword] MULTI}.  No other data may be provided on this line,
though a comment, preceded by the comment marker, is valid.  A warning shall
be produced when a key which has not been specified to be non-unique is
repeated; in this case, the former value shall be overwritten if
\code{overwrite} is \code{true}, otherwise the line shall be ignored and
counted as one that could not be parsed.  It should be noted that non-unique
keys may be of mixed type (even the \code{TYPE} and \code{METADATA} complex
types). For example:
\begin{verbatim}
comment     MULTI   # a comment
comment     STR     some string
comment     F32     1.23456
comment     BOOL    T
\end{verbatim}

If a line does not conform to the rules laid out here, a warning shall
be generated, it shall be ignored and counted as a line that could not
be parsed.  The total number of lines that were not able to be parsed
(including those that were ignored because \code{overwrite} is
\code{false}, and any other parsing problems, but not including blank
lines and comment lines) shall be returned by the function in the
argument \code{nFail}.

Here are some examples of lines of a valid configuration file:
\filbreak
\begin{verbatim}
Double     F64     1.23456789      # This is a comment
Float    F32 0.98765 # This is a comment too
String  STR This is the string that forms the value #comment

 # This is a comment line and is to be ignored
boolean     BOOL    T # The value of `boolean' is `true'

@primes U8  2,3 5 7,11,13 17 #   These are prime numbers

comment MULTI # The rest of this line is ignored, but `comment' is set to be non-unique
comment STR This
comment STR     is
comment STR       a
comment STR        non-unique
comment STR                  key
Float F64 1.23456 # This generates a warning, and, if `overwrite' is `false', is ignored
\end{verbatim}

Of course, a real configuration file should look much nicer to humans
than the above example, but PSLib must be able to parse such ugly
files.

\paragraph{Complex Types}
\subparagraph{TYPE}

We support a modest tree structure by defining a reserved keyword \code{TYPE}.
Any line in the config file which starts with the word \code{TYPE} shall be
interpreted as defining a new valid type.  The defined type name follows the
word \code{TYPE}, and is in turn followed by an arbitrary number of words.
These words are to be interpreted as the names of an embedded \code{psMetadata}
entry, where the values are given on any line which (following the \code{TYPE}
definition) employs the new type name.  For example, a new type may be defined
as:
\begin{verbatim}
TYPE      CELL   EXTNAME   BIASSEC  CHIP
CELL.00   CELL   CCD00     BSEC-00  CHIP.00
CELL.01   CELL   CCD01     BSEC-01  CHIP.00
\end{verbatim}

When \code{psMetadataConfigParse} encounters the \code{TYPE} line, it
should construct a \code{psMetadata} container and fill it with
\code{psMetadataItems} having the names \code{EXTNAME, BIASSEC, CHIP},
with type \code{PS_DATA_STR}, but data allocated.  When it next
encounters an entry of type \code{CELL}, it should then use the given
name (e.g., \code{CELL.00}) for the \code{psMetadataItem}, and copy
the \code{psMetadata} data onto the \code{psMetadataItem.data.md}
entry, filling in the values from the rest of the line (\code{CCD00,
BSEC-00, CHIP.00}).  This hierarchical structure is illustrated in
Figure~\ref{fig:metadata}.

\subparagraph{METADATA}

Another way to form a tree-like structure is to directly define a
\code{psMetadata} entry using a sequence of successive lines to define the
values of the \code{psMetadataItem} entries.  The initial line defines the new
\code{psMetadata} entry and its name.  The following lines have the same format
as the other metadata config file entries.  The sequence is terminated with a
line with a single word \code{END}.  For example, a metadata entry may be
defined as:
\begin{verbatim}
CELL      METADATA
 EXTNAME   STR   CCD00
 BIASSEC   STR   BSEC-00
 CHIP      STR   CHIP.00
 NCELL     S32   24
END
\end{verbatim}

\paragraph{Scoping Rules}

A simple set of ``Scoping Rules'' are required to properly parse a
configuration file.  ``Scope'' refers to the current ``level'' of
\code{METADATA} that a statement appears in.  Statements that are not contained
in a nested \code{METADATA} are said to be in the ``Top level scope''.  Each
level of nested \code{METADATA} statements create a new ``lower level scope''.

\begin{itemize}
\item 
Variable names are unique only to the current level of scope.

\item
non-unique keywords (\code{MULTI}) apply only to the current scope.  i.e. They
are invalid in ``higher'' or ``lower'' level scopes.

\item
\code{TYPE} declarations apply only to the current scope.

\item
\code{METADATA} declarations must begin and end in the same scope.  i.e.  They
may not be declared and end in two different nested METADATA and the same
depth.
\end{itemize}

A series of test inputs is contained in
\S\ref{sec:configtest}.

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

\subsection{Lookup Tables}

Lookup tables store a variety of values indexed on a certain column.
An example is for storing the difference between UT1 and UTC, and the
polar motion vector as a function of date.

One of the key functionalities of a lookup table is to read data from
an ordinary text file into an array of vectors.  This functionality is
generally useful, and so we specify a separate function that may be
called independently:
\begin{prototype}
psArray *psVectorsReadFromFile(const char *filename, const char *format);
\end{prototype}
\code{psVectorsReadFromFile} shall return an array of
\code{psVector}s, read from the specified \code{filename}.  The file
shall be plain text, consisting of an identical number of columns on
each line, with the values separated by whitespace.  Lines commencing
with a comment character (the pound sign, \code{#}) and blank lines
shall be ignored.  The \code{format} is a \code{scanf}-like format
which specifies the number of columns in the file, as well as their
types.  The following formats shall be defined: \code{\%d} for psS32,
\code{\%ld} for psS64, \code{\%f} for psF32, and \code{\%lf} for
psF64.  A star (\code{*}) in the format shall indicate that the column
is to be skipped.

\begin{datatype}
typedef struct {
    const char *filename;               ///< File from which data is to be read
    const char *format;                 ///< scanf-like format string for file
    long indexCol;                      ///< Column of the index vector (starting at zero)
    psVector *index;                    ///< Index values
    psArray *values;                    ///< Corresponding values: an array of vectors
    const double validFrom;             ///< Minimum index value for validity
    const double validTo;               ///< Maximum index value for validity
} psLookupTable;
\end{datatype}

\code{filename} shall specify the file from which the lookup table
data is to be read.  \code{format} shall contain a \code{scanf}-like
format string specifying how the columns are to be interpreted (see
\code{psVectorsReadFromFile}).  \code{indexCol} shall specify the
index of the column (with the first column having an index of zero)
that will form the index values.  \code{index} shall contain the index
values, which shall be sorted in increasing order.  The \code{values}
shall consist of an array of vectors, each of the same length as the
\code{index} vector.  The vectors (including the \code{index} and all
vectors in the \code{values} array) may be any numerical type except
complex types.  The \code{validFrom} and \code{validTo} shall specify
the range of valid values for the index; in most cases, these will
simply be the first and last indices.

The constructor shall be:
\begin{prototype}
psLookupTable *psLookupTableAlloc(const char *filename, ///< File from which to read
                                  const char *format,   ///< scanf-like format string
                                  long indexCol         ///< Column of the index vector (starting at zero)
                                  );
\end{prototype}
This function shall allocate a \code{psLookupTable}, and set the
appropriate values, but it shall not read the lookup table.  This is
so that the lookup table can be specified at the initialisation of a
program, but not read unless required.

The destructor shall free all the components.

\begin{prototype}
psLookupTable *psLookupTableImport(psLookupTable *table,    ///< Lookup table into which to import
                                   const psArray *vectors,  ///< Array of vectors
                                   long indexCol            ///< Index of the index vector in the array of vectors
                                   );
\end{prototype}
\code{psLookupTableImport} shall import an array of vectors into a
\code{table}.  If \code{table} is \code{NULL}, a new
\code{psLookupTable} shall be allocated and returned.  The array of
\code{vectors}, which was likely generated by
\code{psVectorsReadFromFile}, are imported by setting the
\code{table->index} to the vector specified by the \code{indexCol},
and pointing the \code{table->values} array data to the remaining
vectors in \code{vectors}.  Reference counters for the vectors shall
be incremented as appropriate.  The \code{validFrom} and
\code{validTo} members of the \code{table} shall be set to the first
and last values in the index vector.  If the \code{index} vector is
not sorted in the file, the lookup table shall be sorted prior to the
function returning.

\begin{prototype}
long psLookupTableRead(psLookupTable *table);
\end{prototype}
\code{psLookupTableRead} combines \code{psVectorsReadFromFile} and
\code{psLookupTableImport} to read the appropriate file and import the
data into the extant \code{table}.  If the input \code{table} has
already been read from a file, the file shall be re-read, and the
contents replaced.  The function shall return the number of lines read
(not including ignored lines).

Interpolation on a lookup table is performed by the following
functions:

\begin{prototype}
double psLookupTableInterpolate(const psLookupTable *table, double index,
                                long column);
psVector *psLookupTableInterpolateAll(const psLookupTable *table,
                                      double index, psVector *stats);
\end{prototype}
Both functions shall interpolate the \code{table} at the provided
\code{index}.  For \code{psLookupTableInterpolate}, only the value in
the specified \code{column} shall be calculated and returned.  For
\code{psLookupTableInterpolateAll}, all the values shall be calculated
and returned as a \code{psVector}, the type of which shall be
\code{PS_TYPE_F64}.

\tbd{What's the stats for???}

If the \code{index} is beyond the range of the \code{table},
\code{psLookupTableInterpolate} shall return \code{NaN}, and
\code{psLookupTableInterpolateAll} shall return \code{NULL} --- that
is, no attempt is made at extrapolation.

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

\section{Mathematical Structures}
\label{sec:mathStructures}

Throughout PSLib, we require a variety of structures which correspond
to different mathematical data concepts.  For example, we have a data
structure which corresponds to one-dimensional arrays (vectors) of
different data types (\code{int}, \code{float}, etc).  We also have a
data structure which corresponds to two-dimensional arrays (images or
matrices), again with different data types for the individual
elements.

A variety of functions perform operations which are equivalent for
different data types of the same dimension, or may even be defined for
different data types of different dimensions.  For example, if we
write the operation $x + y$, the operation is clearly defined
regardless of whether the operands $x$ and $y$ are both zero
dimensional (single numbers), one dimensional (vectors), two
dimensional (images), etc. It is even reasonable to define the meaning
of such an operation if the data dimensions do not match: if $x$ is a
scalar and $y$ is an image, the natural operation is to add the value
of $x$ to every element of $y$; we can also define the meaning of the
operation if $x$ is a vector and $y$ is a matrix.  Nor does it matter
mathematically that the element data types match; the sum of a float
and an integer is a well-defined quantity.  One constraint should be
noted: the size of the elements in each dimension must match.  For
example, if $x$ were a vector of 100 elements, but $y$ were a vector
of 1000 elements, the meaning of the operation $x + y$ is unclear.
This type of operation should be invalid and should generate an error.

Given that some functions should be able to operate equivalently (or
identically) on a wide range of data types, we define a mechanism
which allows the C functions to accept a generic data type, and
determine the type of the data on the basis of the data.  
Supported data types must be defined by a structure in which
the first element is always of type \code{psMathType}:
\begin{datatype}
typedef struct {
    psDimen dimen;                      ///< The dimensionality
    psElemType type;                    ///< The type
} psMathType;
\end{datatype}
where \code{psDimen dimen} defines the dimensionality of the data:
\begin{datatype}
typedef enum {
    PS_DIMEN_SCALAR,                    ///< Scalar
    PS_DIMEN_VECTOR,                    ///< A vector
    PS_DIMEN_TRANSV,                    ///< A transposed vector
    PS_DIMEN_IMAGE,                     ///< An image (matrix)
    PS_DIMEN_OTHER                      ///< Not supported for arithmetic
} psDimen;
\end{datatype}
\code{psElemType type}, which defines the data type of each element, has
already been defined

We discuss the application of \code{psMathType} in more detail in
section~\ref{sec:arithmetic}.  

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

\subsection{Math Casting}

We define a basic data type which only contains the type information.
This structure should be used to cast an unknown \code{psMath} data
structure (\code{psImage}, \code{psVector}, \code{psScalar}) so the
data type testing may be safely performed.  

\begin{datatype}
typedef struct { 
    psMathType type;                    ///< data type information 
} psMath;
\end{datatype}

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

\subsection{Scalars}

We define a basic scalar data type which includes the type
information.  This structure allows scalars to be used in functions
which interpret the data type from the structure when deciding how to
perform an operation.  The basic scalar structure is:
\begin{datatype}
typedef struct {
    psMathType type;                    ///< data type information
    union {                            
        psS8   S8;                      ///< bye value entry
        psS16 S16;                      ///< short int value entry
	psS32 S32;                      ///< int value entry
	psS64 S64;                      ///< long int value entry
        psU8   U8;                      ///< unsigned byte value entry
        psU16 U16;                      ///< unsigned short int value entry
	psU32 U32;                      ///< unsigned int value entry
	psU64 U64;                      ///< unsigned long int value entry
        psF32 F32;                      ///< float value entry
        psF64 F64;                      ///< double value entry
        psC32 C32;                      ///< complex value entry
        psC64 C64;                      ///< double complex value entry
    } data;
} psScalar;
\end{datatype}

In addition, we specify two functions for working with \code{psScalar} data:
\begin{prototype}
psScalar *psScalarAlloc(complex value, psElemType type);
psScalar *psScalarCopy(const psScalar *value);
\end{prototype}
The first creates a \code{psMathType}-ed structure from a constant
value, casting it as appropriate based on the \code{type}.  The second
copies the provided \code{psScalar} value.  This latter function is
necessary to keep a copy of an existing \code{psScalar} value, since
\code{psBinaryOp} and \code{psUnaryOp} are required to free incoming
\code{psScalar} data (see \S\ref{sec:arithmetic}).

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

\subsection{Vectors}

We require several related types of basic one-dimensional arrays:
arrays of values of type \code{int}, \code{float}, \code{double},
\code{complex float}, and \code{complex}.  We have defined a
single structure, \code{psVector} to represent these concepts:
%
\begin{datatype}
typedef struct {
    psMathType type;                    ///< vector data type and dimension
    long n;                             ///< size of vector
    const long nalloc;                  ///< allocated data block
    union {
        psS8  *S8;                      ///< Pointers to byte data
        psS16 *S16;                     ///< Pointers to short-integer data
        psS32 *S32;                     ///< Pointers to integer data
        psS64 *S64;                     ///< Pointers to long-integer data
        psU8  *U8;                      ///< Pointers to unsigned-byte data
        psU16 *U16;                     ///< Pointers to unsigned-short-integer data
        psU32 *U32;                     ///< Pointers to unsigned-integer data
        psU64 *U64;                     ///< Pointers to unsigned-long-integer data
        psF32 *F32;                     ///< Pointers to floating-point data
        psF64 *F64;                     ///< Pointers to double-precision data
        psC32 *C32;                     ///< Pointers to complex floating-point data
        psC64 *C64;                     ///< Pointers to complex floating-point data
    } data;
    void *lock;                         ///< Lock for thread safety
} psVector;
\end{datatype}
%
In this structure, the argument \code{n} is the length of the array
(the number of elements); \code{nalloc} is the number of elements
allocated (\code{nalloc} $\ge$ \code{n}).  The allocated memory is
available in the union \code{data} which consists of pointers to each
of the defined primitive data types.  Note the parallelism in the
names of the types, union elements, and the \code{psElemType} names.
This parallelism allows us to use automatic construction mechanisms
effectively.  The data type is defined by the first element,
\code{psMathType}.  The structure is associated with a constructor and
reallocator:
%
\begin{prototype}
psVector *psVectorAlloc(long nalloc, psElemType type);
psVector *psVectorRealloc(psVector *vector, long nalloc);
psVector *psVectorRecycle(psVector *vector, long nalloc, psElemType type);
\end{prototype}
%
In these functions, \code{nalloc} is the number of elements to
allocate.  For \code{psVectorAlloc}, the value of \code{psVector.n} is
set to \code{nalloc}.  Users may choose to restrict the data range
after the \code{psVectorAlloc} function is called.  For
\code{psVectorRealloc}, if the value of \code{nalloc} is smaller than
the current value of \code{psVector.n}, then \code{psVector.n} is set
to \code{nalloc}, the array is adjusted down to match \code{nalloc},
and the extra elements are lost.  If \code{nalloc} is larger than the
current value of \code{psVector.n}, \code{psVector.n} is left intact.
If the value of \code{vector} is \code{NULL}, then
\code{psVectorRealloc} must generate an error.

\code{psVectorRecycle} shall recycle the input \code{vector}, such
that the output \code{psVector} matches the length required for
\code{nalloc} elements of the specified \code{type}.  In the event
that the input \code{vector} is \code{NULL}, a new \code{psVector}
shall be allocated and returned.

\begin{prototype}
psVector *psVectorExtend(psVector *vector, long delta, long nExtend);
\end{prototype}

This function increments \code{psVector.n} (the number of elements in
the vector) by \code{nExtend}.  If the current length of the vector
plus {\em twice} the number of new elements (\code{nExtend}) is
greater than the allocated space, an additional \code{delta} elements
are allocated.  If the value of \code{delta} is less than 1, 10 shall
be used.

Here is an example of how \code{psVectorExtend} is used to
automatically increment the vector length.
\begin{verbatim}
  // create data vector
  psVector *y = psVectorAlloc(100);
  y->n = 0;
  for (int i = 0; i < 1000; i++) {
    y->data.F32[y->n + 0] = 2*i;
    y->data.F32[y->n + 1] = 2*i;
    y->data.F32[y->n + 2] = 2*i;
    y = psVectorExtend (y, 100, 3);
    // increments n by 3, extends length by 100 if needed
  }
\end{verbatim}
Note that the specification that the allocation always be greater than
the number of elements by twice the number of new elements implies
that there will be room on the next loop for \code{nExtend} new
elements, as in this example.

\begin{prototype}
psVector *psVectorCopy(psVector *output, const psVector *input, psElemType type);
\end{prototype}

\code{psVectorCopy} shall copy the elements in the \code{input} vector
to the \code{output} vector, casting to the specified \code{type}.  In
the event that the \code{output} is \code{NULL}, a new \code{psVector}
shall be allocated.  The returned \code{psVector} shall be of the
given \code{type}.

\begin{prototype}
bool psVectorInit(psVector *vector, ...);
\end{prototype}

\code{psVectorInit} shall initialize the vector with the given value.
The input data is cast to match the vector datatype, allowing for
integers to be preserved.

\begin{prototype}
bool psVectorSet(const psVector *input, long position, complex value);
complex psVectorGet(const psVector *input, long position);
\end{prototype}

These accessor functions are provided as a convenience to the user.
\code{psVectorSet} sets the value of the \code{input} vector at the
specified \code{position} to \code{value} (appropriately cast),
returning \code{true} if successful.  \code{psVectorGet} returns the
value of the \code{input} vector at the specified \code{position}.  A
negative \code{position} means index from the end.

\begin{prototype}
psVector psVectorCreate(psVector *input, double lower, double upper, double delta, psElemType type);
\end{prototype}

This function creates a new vector, or reallocates the provided vector
if \code{input} is not \code{NULL}.  The created vector consists of
the data range starting at \code{lower}, running to \code{upper}, in
steps of \code{delta}.  The upper-end value is {\em exclusive}; the
sequence is equivalent to \code{for (x = lower; x < upper; x +=
  delta)}.  

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

\subsection{Images}

The most important data product produced by the telescope is an image.
The simplest image is a 2-D collection of pixels, each with some
value.  We require a basic image data type:

\begin{datatype}
typedef struct psImage {
    psMathType type;                    ///< image data type and dimension
    const int numCols;                  ///< Number of columns in image
    const int numRows;                  ///< Number of rows in image.
    const int col0;                     ///< Column position relative to parent.
    const int row0;                     ///< Row position relative to parent.
    union {
        psS8  **S8;                     ///< Pointers to char data
        psS16 **S16;                    ///< Pointers to short-integer data
        psS32 **S32;                    ///< Pointers to integer data
        psS64 **S64;                    ///< Pointers to long-integer data
        psU8  **U8;                     ///< Pointers to unsigned-char data
        psU16 **U16;                    ///< Pointers to unsigned-short-integer data
        psU32 **U32;                    ///< Pointers to unsigned-integer data
        psU64 **U64;                    ///< Pointers to unsigned-long-integer data
        psF32 **F32;                    ///< Pointers to floating-point data
        psF64 **F64;                    ///< Pointers to double-precision data
        psC32 **C32;                    ///< Pointers to complex floating-point data
        psC64 **C64;                    ///< Pointers to complex floating-point data
        psPtr V;                        ///< Pointers to untyped data
    } data;
    const struct psImage *parent;       ///< parent, if a subimage 
    psPtr p_rawDataBuffer;              ///< Raw data; private
    psArray *children;                  ///< children of this region
    void *lock;                         ///< Lock for thread safety
} psImage;
\end{datatype}

This structure represents an image consisting of a 2-D array of
pixels.  The size of this array is given by the elements
\code{(numRows, numCols)}.  The data type of the pixel is defined in
the \code{psMathType type} entry (specifically, the \code{psElemType}
member, \code{type}; see \S\ref{sec:arithmetic}).  (N.B.\ that for
FITS images, these values are restricted to the datatypes equivalent
to the valid BITPIX values 8, 16, 32, -32, -64).  The image
represented in the data structure may represent a subset of the pixels
in a complete array, in which case the image is considered to be the
child of that parent array.  The offset of the \code{(0,0)} pixel in
this array relative to the parent array is given by the elements
\code{(col0,row0)}: \code{col0} is the starting column number in the
parent image while \code{row0} is the starting row number.  The
structure may include references to subrasters (\code{children})
and/or to a containing array (\code{parent}).  Unless this image is a
child of another image (represents a subset of the pixels of another
image), the image data is allocated in a contiguous block
(\code{data.V}).

\subsubsection{Support Functions}

We define the following supporting functions for images, which are
valid for data types \code{psS8, psS16, psU8, psU16, psF32, psF64,
psC32, psC64}.

\begin{prototype}
psImage *psImageAlloc(int numCols, int numRows, psElemType type);
\end{prototype}
Create an image of a specified \code{numCols}, \code{numRows}, and data
\code{type}.  This function must allow any of the valid image data
types and not restrict to the valid FITS BITPIX types.  The image
dimensionality must be 2.  

\begin{prototype}
psImage* psImageRecycle(
    psImage* old,                       ///< the psImage to recycle by resizing image buffer
    int numCols,                        ///< the desired number of columns in image
    int numRows,                        ///< the desired number of rows in image
    const psElemType type               ///< the desired datatype of the image
);
\end{prototype}
\code{psImageRecycle} shall recycle the input \code{old} image, such
that the output \code{psImage} matches the specified size
(\code{numCols}$\times$\code{numRows}) and \code{type}.  In the event
that the input \code{old} image is \code{NULL}, a new \code{psImage}
shall be allocated and returned.

\begin{prototype}
int psImageFreeChildren(psImage *image);
\end{prototype}

\code{psImageFreeChildren} shall free all child images of the given
\code{image}.

\begin{prototype}
bool psImageInit(psImage *image, ...);
\end{prototype}

\code{psImageInit} shall initialize the image with the given value.
The input data is cast to match the image datatype, allowing for
integers to be preserved.

\begin{prototype}
bool psImageSet(const psImage *image, int x, int y, complex value);
complex psImageGet(const psImage *image, int x, int y);
\end{prototype}

These accessor functions are provided as a convenience to the user.
\code{psImageSet} sets the value of the \code{image} at the specified
\code{x,y} position to \code{value} (appropriately cast), returning
\code{true} if successful.  \code{psImageGet} returns the value of the
\code{image} at the specified \code{x,y} position.  A negative value
for the \code{x,y} position means index from the end.

\subsubsection{Image Regions}

In many places, we need to refer to a rectangular area.  We define a
structure to represent a rectangle:
\begin{datatype}
typedef struct {
    float x0;
    float x1;
    float y0;
    float y1;
} psRegion;
\end{datatype}
This structure is used in psLib as an abbreviation for the four
floats, and is defined statically.  Functions which accept or return a
\code{psRegion} shall do so by value, not by pointer.

All functions which use a \code{psRegion} must interpret the
definition of $(x0,y0)$ and $(x1,y1)$ in the same way.  The coordinate
$(x0,y0)$ defines the starting pixel in the region.  The coordinate
$(x1,y1)$ defines the outer bound of the region, and are NOT included
in the region.  The size of the region is $(x1-x0,y1-y1)$.  If either
$x1$ or $y1$ is less than or equal to 0, the value is added to the
image dimensions (e.g., $Nx + x1$).  Thus a region \code{[0:0,0:0]}
refers to the full image array, while \code{[0:-10,0:-20]} refers to
the entire image, minus the last 10 columns and the last 20 rows.

We define two functions to set and return the value of a
\code{psRegion}.  The first defines the region by the corner
coordinates.  The second function converts the IRAF description region
in the form \code{[x0:x1,y0:y1]}, used for header entries such as
\code{BIASSEC}, into the corresponding \code{psRegion} structure (any
values that do not parse correctly shall be returned as \code{NaN}).
We also define a function that converts a \code{psRegion} to the
corresponding IRAF description.

\begin{prototype}
psRegion psRegionSet(float x0, float x1, float y0, float y1);
psRegion psRegionFromString(const char *region);
psString psRegionToString(const psRegion region);
\end{prototype}

{\bf Note that regions specified by strings are in the FITS standard.}
It is the responsibility of \code{psRegionFromString} and
\code{psRegionToString} to convert between the PS standard (0 means
first pixel; upper value is not included, but lower is) and the FITS
standard (1 means first pixel; lower and upper values are included),
which simply involves subtracting one from \code{x0} and \code{y0}
when going from a string representation to a \code{psRegion}.

\begin{prototype}
psRegion psRegionForImage(psImage *image, psRegion in);
\end{prototype}

An image region defined with negative upper limits may be rationalized
for the bounds of a specific image with \code{psRegionForImage}.  The
output of this function is a region with negative upper limits
replaced by their corrected value appropriate to the given image.  In
addition, the lower and upper limits are forced to lie within the
bounds of the image.  If the lower-limit coordinates are less than the
lower bound of the image, they are limited to the lower bound of the
image.  Conversely, if the upper-limit coordinates are greater than
the upper bound of the image, they are truncated to define only valid
pixels.  If the lower-limit coordinates are greater than the upper
bounds of the image, or the upper-limit coorindates are less than the
lower bounds of the image, the coordinates should saturate on those
limits.  The output of this function is always a valid region, though
it may define an area of 0 pixels.


\begin{prototype}
psRegion psRegionForSquare(double x, double y, double radius);
\end{prototype}
This utility function defines a \code{psRegion} corresponding to the
square with center at coordinate \code{x,y} and with code{radius}.
The width of the square is thus \code{2radius + 1}.

\pagebreak 

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

\section{Input/Output}

\subsection{XML Functions}

Within Pan-STARRS, we will use XML documents as a transport mechanism
to carry data between programs and between IPP and other subsystems.
Configuration information may be stored as well as XML documents, in
addition to the text format discussed in the discussion on Metadata.
XML is an extremely variable document format, and it is not currently
the intention of PSLib to provide a complete PSLib version of XML
operations.  Rather, a limited number of operations are defined to
convert specific data structures to an appropriate XML document.  To
maximize the simplicity of the XML APIs, we will use the convention
that a single XML document to be parsed by PSLib shall contain only a
single data structure.  Each of the XML APIs therefore take as input a
reference to a complete XML document and return a PSLib data
structure, or take a PSLib data structure and return a complete XML
document.

We start by defining a PSLib wrapper type which is a pointer to an XML
document in memory.  We wrap the \code{libxml2} version of an XML
document pointer for now:
\begin{datatype}
typedef xmlDocPtr psXMLDoc;
\end{datatype}

\begin{prototype}
psXMLDoc *psXMLDocAlloc(void);
\end{prototype}

The next pair of functions convert a \code{psMetadata} data structure
to a complete \code{psXMLDoc} (in memory) and vice versa:
\begin{prototype}
psXMLDoc *psMetadataToXMLDoc(const psMetadata *md);
psMetadata *psXMLDocToMetadata(const psXMLDoc *doc);
\end{prototype}

The next pair of functions loads the data in a named file into a
complete \code{psXMLDoc} (in memory) and write out the \code{psXMLDoc}
to a named file:
\begin{prototype}
psXMLDoc *psXMLParseFile(const char *filename);
bool psXMLDocToFile(const psXMLDoc *doc, const char *filename);
\end{prototype}

The next pair of functions accepts a block of memory and parses it
into a complete \code{psXMLDoc} (also in memory), and vice versa:
\begin{prototype}
psXMLDoc *psXMLParseMem(const char *buffer, int size);
bool psXMLDocToMem(const psXMLDoc *doc, char *buffer);
\end{prototype}

The next pair of functions read from and write to a file descriptor.
The first converts the incoming data to a complete \code{psXMLDoc}
(also in memory), the second writes the \code{psXMLDoc} to the file
descriptor:
\begin{prototype}
psXMLDoc *psXMLParseFD(int fd);
bool psXMLDocToFD(const psXMLDoc *doc, int fd);
\end{prototype}

\subsection{Database Functions}

Many of the applications that PSLib will be used for will require
access to a simple relational database.  PSLib includes generic
database-independent interface mechanisms as part of its API set.  The
most important aspect of PSLib's database support is to abstract as
much database specific complexity as is feasible.  As almost all RDBMS
provide at least a simple transactional model, commit and rollback
support should be provided.

Currently, only support for MySQL 4.1.x is required but other backends
may be added as options in the future.  As a particular example which
has implications for the database interaction model, support for
SQLite may be required in the future.  Currently, the choice of
backend database interface may be made as a compile option.  Details
of the specified APIs in the discussion below refer to the relevant
MySQL functions.

Database errors must be trapped and placed onto the psError stack.
The complete error message should be retrieved with the database's
error function.

\subsubsection{Managing the Database Connection}

We specify a database handle which carries the information about the
database connection:

\begin{datatype}
typedef struct {
    void *mysql;
} psDB;
\end{datatype}

The following collection of functions provides basic database functionality:

\begin{prototype}
psDB *psDBInit(const char *host, const char *user, const char *passwd, const char *dbname);
void psDBCleanup(psDB *dbh);
bool psDBCreate(psDB *dbh, const char *dbname);
bool psDBChange(psDB *dbh, const char *dbname);
bool psDBDrop(psDB *dbh, const char *dbname);
\end{prototype}

For MySQL support, \code{psDBInit} wraps \code{mysql_init} and
\code{mysql_real_connect} in order to initialize a psDB structure and
establish a database connection.  A null pointer should be returned on
failure.

When implementing support for SQLite, or other DB which is purely
file-based, the \code{host}, \code{user}, and \code{passwd} arguments
would be ignored while \code{dbname} would specify the path to the
SQLite db file.

\code{psDBCleanup} shall wrap \code{mysql_close}.  \code{psDBCreate}
shall wrap \code{mysql_create_db}.  \code{psDBChange} shall wrap
\code{mysql_select_db}.  \code{psDBDrop} shall wrap
\code{mysql_drop_db}.

\subsubsection{Interacting with Database Tables}

The functions in this section perform high level interactions with the
database tables.  All of them should behave ``atomically'' with
respect to the state of the database.  Specifically, all interactions
with the database should be done as a part of a transaction that is
rolled-back on failure and committed only after all queries used by
the API have been run.  In general, this API set attempts to treat a
database table as a 2D matrix where columns can be represented by a
\code{psVector} and rows as a \code{psMetadata} type.  A
\code{psMetadata} collection is also used to define the columns of a
table and as part of the query restrictions.

\begin{prototype}
bool psDBCreateTable(psDB *dbh, const char *tableName, const psMetadata *md);
\end{prototype}

This function generates and executes the SQL needed to create a table
named \code{tableName}, with the column names and datatypes as
described in \code{md}.  Each data item in the \code{psMetadata}
collection represents a single table field.  The name of the field is
given by the name of the \code{psMetadataItem} and the data type is
given by the \code{psMetadataItem.type} entry.  A lookup table should
be used to convert from PSLib types into MySQL compatible SQL data
types.  For example, a \code{PS_DATA_STR} would map to an SQL99
varchar.  If the value of \code{type} is \code{PS_DATA_STR} then the
\code{psMetadataItem.data} element is set to a string with the length
for the field written as a text string.  The value of the
\code{psMetadataItem.data} element is unused for the
\code{PS_DATA_PRIMITIVE} types.  Other metadata types beyond
\code{PS_DATA_STR} and \code{PS_DATA_PRIMITIVE} are not allowed in a
table definition metadata collection.

Database indexes can be specified by setting the \code{comment} field to
``\code{Primary Key}'' or ``\code{Key}''.  ``Auto-incrementing'' columns may be
specified by setting the \code{comment} field to ``\code{AUTO_INCREMENT}''.
Indexes and auto-increments maybe be combined in the same \code{comment}.  They
must be separated by optional whitespace, a comma, and then more optional
whitespace.  The \code{comment} field is otherwise ignored.

\begin{prototype}
bool psDBDropTable(psDB *dbh, const char *tableName);
\end{prototype}

This function deletes the specified table.

\begin{prototype}
bool p_psDBRunQuery(psDB *dbh, const char *format, ...);
\end{prototype}

This function will execute a string as a raw SQL query.  \code{format} is a
\code{printf} style formatting code to be implimented with \code{vsprintf()}.
No additional processing of the string or abstraction of the underlying
database's SQL dialect is provided.

\begin{prototype}
psArray *psDBSelectColumn(psDB *dbh, const char *tableName, const char *col, unsigned long long limit);
psVector *psDBSelectColumnNum(psDB *dbh, const char *tableName, const char *col,
                              psElemType type, unsigned long long limit);
\end{prototype}

These functions generates and executes the SQL needed to select an entire
column from a table or up to \code{limit} rows from it.  If \code{limit} is 0,
the entire range is returned.  The database response is processed and a
\code{psArray} of strings is returned.  The Num version of the function returns
the data in a \code{psVector}, data cast to \code{type}.  It returns an error
(NULL) if the requested field is not a numerical type.

\begin{prototype}
psArray *psDBSelectRows(psDB *dbh, const char *tableName, const psMetadata *where, unsigned long long limit);
\end{prototype}

This function returns rows from the specified table which match
the restrictions given by \code{where}.  The restrictions are
specified as field / value pairs.  The \code{psMetadata} collection
where must consist of valid database fields, though the database query
checking functions may be used to validate the fields as part of the
query.  If \code{where} is \code{NULL}, then there are no restrictions
on the rows selected.  The selected rows are returned as a
\code{psArray} of \code{psMetadata} values, one per row. 

\begin{prototype}
bool psDBInsertOneRow(psDB *dbh, const char *tableName, const psMetadata *row);
\end{prototype}

Insert the data from \code{row} into \code{tableName}.  It should be noted in
the API reference that if fields are specified in \code{row} that do not exist
in \code{tablename}, the insert will fail.

\begin{prototype}
bool psDBInsertRows(psDB *dbh, const char *tableName, const psArray *rowSet);
\end{prototype}

Similar to \code{psDBInsertOneRow()}, this function inserts many rows at once
and is atomic for the complete set of rows.

\begin{prototype}
psArray *psDBDumpRows(psDB *dbh, const char *tableName);
\end{prototype}

Fetch all rows as an psArray of psMetadata.

\begin{prototype}
psMetadata *psDBDumpCols(psDB *dbh, const char *tableName);
\end{prototype}

Fetch all columns, as either a psVector or a psArray depending on whether or not
the column is numeric, and return them in a psMetadata structure where
psMetadataItem.name contains the column's name.

\begin{prototype}
long psDBUpdateRows(psDB *dbh, const char *tableName, const psMetadata *where,
                    const psMetadata *values);
\end{prototype}

Update the columns contained in \code{values} in the row(s) that have a field
with the value indicated by \code{where} (note that this is only allows very
limited use of SQL99's ``where'' semantics).  The number of rows modified is
returned.  A negative value is return to indicate an error. If there are
multiple psMetadataItems in \code{where} then each item should be considered as
an additional constraint.  e.g.  ``where foo = x and where bar = y''

\begin{prototype}
long psDBDeleteRows(psDB *dbh, const char *tableName, const psMetadata *where, unsigned long long limit);
\end{prototype}

Delete the rows that are matched by \code{where} using the same semantics for
\code{where} as in psDBUpdateRow().  \code{limit} specifies the maximum number
of rows that may be deleted per invocation.  A negative value is returned to
indicate an error.

\subsection{FITS I/O Functions}

We need a variety of I/O functions between the disk and certain of our
PSLib data structures.  We need the ability to access FITS headers,
images and tables (both ASCII and Binary).  We define here the FITS
I/O functions, all of which are currently specified as wrappers to
functions within CFITSIO.  CFITSIO provides a wide range of utilities
which we do not feel are particularly appropriate as part of a generic
I/O library, such as assumptions about names which change the data
interpretation, etc.  We are defining our calls to avoid the hidden
'features'.  The CFITSIO functions which are wrapped should in general
be the most basic versions.

\begin{datatype}
typedef struct {
    fitsfile fd;                   // cfitsio structure
    bool writable;                 // Is the file writable?
} psFits;
\end{datatype}
We begin by defining a datatype to wrap the CFITSIO \code{fitsfile}
structure.  This is necessary to allow repeated access to the data in
a file without multiple open commands (which are expensive).  Write
operations are only permitted if \code{writable} is \code{true}.

\subsubsection{FITS File Manipulations}

\begin{prototype}
psFits *psFitsOpen(const char *filename, const char mode);
\end{prototype}

Opens a FITS file and positions the pointer to the PHU.  The file is
opened in the requested \code{mode}, which may be one of \code{r}
(read only) \code{r+} (read and write), \code{rw} (alias for
\code{r+}) or \code{w} (create new file for writing).

\begin{prototype}
psFits *psFitsOpenFD(int fd);
\end{prototype}

Opens a file descriptor and positions the pointer to the PHU.

\begin{prototype}
psFits *psFitsOpenStream(FILE *stream);
\end{prototype}

Opens a stream and positions the pointer to the PHU.

\begin{prototype}
bool psFitsClose(psFits *fits);
\end{prototype}

Closes a FITS file.

\begin{prototype}
bool psFitsMoveExtName(const psFits *fits, const char *extname);
\end{prototype}

Positions the pointer to the beginning of the specified
\code{extname}.  If the \code{extname} does not exist, the function
shall fail.  

\begin{prototype}
bool psFitsMoveExtNum(const psFits* fits, int extnum, bool relative);
\end{prototype}

Moves the pointer to the beginning of the specified HDU number.  If
\code{relative} is TRUE, \code{extnum} represents the number of HDUs
relative to the current HDU.  The PHU is entry number 0, while the
extended data segments start at number 1.

\begin{prototype}
bool psFitsMoveEnd(psFits *fits);
\end{prototype}

Moves the pointer to the end of the file, so that a new extension may
be written.

\begin{prototype}
int psFitsGetExtNum(const psFits* fits);
\end{prototype}

Returns the current HDU number (i.e., file position).  

\begin{prototype}
int psFitsGetSize(const psFits* fits);
\end{prototype}

Returns the number of HDUs in the file.

\begin{prototype}
bool psFitsDeleteExtNum(psFits *fits, int extnum, bool relative);
bool psFitsDeleteExtName(psFits *fits, const char *extname);
\end{prototype}

These functions are similar to \code{psFitsMoveExtNum} and
\code{psFitsMoveExtName} except that they delete the specified
extension.  The file pointer is left pointing to where the extension
was before deletion.

\begin{prototype}
bool psFitsTruncate(psFits *fits);
\end{prototype}

Deletes all extensions after the position of the file pointer.

\begin{datatype}
typedef enum {
    PS_FITS_TYPE_NONE,
    PS_FITS_TYPE_IMAGE,
    PS_FITS_TYPE_BINARY_TABLE,
    PS_FITS_TYPE_ASCII_TABLE,
    PS_FITS_TYPE_ANY
} psFitsType;
\end{datatype}

\begin{prototype}
psFitsType psFitsGetExtType(const psFits* fits);
\end{prototype}

Gets the current HDU's type (table or image).

\begin{prototype}
psString psFitsGetExtName(const psFits* fits);
bool psFitsSetExtName(psFits* fits, const char* name);
\end{prototype}

\code{psFitsGetExtName} shall return the name of the current extension
for the given \code{fits} file (as specified by the \code{EXTNAME}
header).  \code{psFitsSetExtName} shall change the name of the current
extension for the given \code{fits} file to \code{name}, returning
\code{true} upon success and \code{false} otherwise.

\subsubsection{FITS Header I/O Functions}

\begin{prototype}
psMetadata *psFitsReadHeader(psMetadata *out, const psFits *fits);
\end{prototype}
Read header data into a \code{psMetadata} structure.  The data is read
from the current HDU pointed at by the \code{psFits *fits} entry.  If
\code{out} is \code{NULL}, a new psMetadata is created.

\begin{prototype}
psMetadata *psFitsReadHeaderSet(const psFits *fits);
\end{prototype}
Load a complete set of headers from the \code{psFits} file pointer.
This function loads the headers from all extensions into a
\code{psMetadata} collection, each entry of which is a pointer to a
\code{psMetadata} structure containing the header data.  The metadata
entry names are the \code{EXTNAME} values for each header (with the
value of \code{PHU} for the primary header unit).  At the start of the
operation, the file pointer is rewound to the beginning of the file.
At the end, it is positioned where it started when the function was
called.

\begin{prototype}
bool psFitsWriteHeader(const psMetadata *output, psFits *fits);
\end{prototype}
Write metadata into the header of a FITS image file.  The header is
written at the current HDU.

\begin{prototype}
psMetadata *psFitsHeaderFromImage(psMetadata *header, psImage *image, const char *extname);
psMetadata *psFitsHeaderFromTable(psMetadata *header, psArray *table, const char *extname);
\end{prototype}

These functions update the given \code{header} to be compatible with
the supplied \code{image} (updating \code{BITPIX, NAXIS1, NAXIS2} or
\code{table} (an array of \code{psMetadata}), and change the
\code{EXTNAME} header to \code{extname}.

\begin{prototype}
bool psFitsHeaderValidate(psMetadata *header);
\end{prototype}

\code{psFitsHeaderValidate} shall validate the supplied \code{header}
so that it is in compliance to the FITS standard for header keyword
names and types.  The contents should be changed to conform to the
FITS standard where possible by truncating names and adding (where
missing) those keywords required by FITS (e.g., \code{SIMPLE}).  If
the resulting \code{header} conforms to the FITS standard, the
function shall return \code{true}; otherwise the function shall return
\code{false}.

\subsubsection{FITS Image I/O Functions}

\begin{prototype}
psImage *psFitsReadImage(psImage *out, const psFits *fits, psRegion region, int z);
\end{prototype}
Read an image or subimage from the \code{psFits} file pointer.  This
function is a wrapper to the CFITSIO library function.  The input
parameters allow a full image or a subimage to be read.  The region to
be read is specified by \code{region}.  A negative value for either of
\code{region.x1} or \code{region.y1} specifies the size of the region
to be read counting down from the end of the array.  

If the native image is a cube, the value of z specifies the requested
slice of the image.  This function must call \code{psError} and return
\code{NULL} if any of the specified parameters are out of range for
the data in the image file, or if the image on disk is zero- or
one-dimensional.  This function need only read images of the native
FITS image types (\code{psU8}, \code{psS16}, \code{psS32},
\code{psF32}, \code{psF64}).  The user is expected to convert the data
type as needed with \code{psImageCopy}.
 
\begin{prototype}
bool psFitsUpdateImage(psFits *fits, const psImage *input, psRegion region, int z);
\end{prototype}
Write an image section to the open \code{psFits} file pointer.  This
operation may write a portion of an image over the existing bytes of
an existing image.  Care must be taken to interpret \code{region},
which specified the output pixels to be written / over-written.  If
the combination of \code{region} and the size of \code{psImage *input}
implies writing pixels outside the existing data area of the image,
the function shall return an error (ie, if \code{region.x0 + image.nx
>= NAXIS1}, \code{region.y0 + image.ny >= NAXIS2}, or \code{z >=
NAXIS3}).  This function will only write images of the native FITS
image types (\code{psU8}, \code{psS16}, \code{psS32}, \code{psF32},
\code{psF64}).  The user is expected to convert the data type as
needed with \code{psImageCopy}.  The return value must be 0 for a
successful operation and 1 for an error.

\begin{prototype}
bool psFitsWriteImage(psFits *fits, psMetadata *header, const psImage *input, int depth,
                      const char *extname);
\end{prototype}
Create a new image based on the dimensions specified for the image and
the requested depth.  The header and image data segments are written
in the file at the current position of the \code{psFits} pointer.
This function will only write images of the native FITS image types
(\code{psU8}, \code{psS16}, \code{psS32}, \code{psF32}, \code{psF64}).
The user is expected to convert the data type as needed with
\code{psImageCopy}.  The return value must be 0 for a successful
operation and 1 for an error.  If \code{extname} is not \code{NULL},
the \code{EXTNAME} in the \code{header} shall be updated to be
\code{extname} before writing.

\subsubsection{FITS Table I/O Functions}

\begin{prototype}
psMetadata *psFitsReadTableRow(const psFits *fits, int row);
\end{prototype}
This function reads a single row of the table in the extension pointed
at by the \code{psFits} file pointer.  The row number to be read is
given by \code{row}.  The result is returned as a \code{psMetadata}
collection with elements of the appropriate types and keys
corresponding to the table column names.  The function must apply the
needed byte-swapping on the data in the row based on the description
of the table data in the table header.  \tbr{we may need to be more
flexible here: if we call this function repeatedly, it would be more
efficient to pass the corresponding header or keep it somewhere (and
the file pointer location, for that matter).}

\begin{prototype}
psPtr psFitsReadTableRowRaw(size_t *size, const psFits *fits, int row);
\end{prototype}
This function reads a single row of the table in the extension pointed
at by the \code{psFits} file pointer.  The row number to be read is
given by \code{row}.  The result is returned as collection of
\code{size} bytes allocated by the function.  The function must
apply the needed byte-swapping on the data in the row based on the
description of the table data in the table header.  \tbr{we may need
to be more flexible here: if we call this function repeatedly, it
would be more efficient to pass the corresponding header or keep it
somewhere (and the file pointer location, for that matter).}

\begin{prototype}
psArray *psFitsReadTableColumn(const psFits *fits, const char *colname);
\end{prototype}
This function reads a single column of the table in the extension
pointed at by the \code{psFits} file pointer.  The column is specified
by the FITS table column key given by \code{row}.  The result is
returned as a \code{psArray}, with the data from one row of the table
column per array element.

\begin{prototype}
psVector *psFitsReadTableColumnNum(const psFits *fits, const char *colname);
\end{prototype}
This function reads a single column of the table in the extension
pointed at by the \code{psFits} file pointer.  The column is specified
by the FITS table column key given by \code{row} and must be of a
numeric data type.  The result is returned as a \code{psVector} of the
appropriate data type, with the data from one row of the table column
per array element.

\begin{prototype}
psArray *psFitsReadTableRaw(size_t *size, const psFits *fits);
\end{prototype}
This function reads the entire data block from a table into the a
\code{psArray}, with one element of the array per row.  The number of
bytes per row is returned in \code{size}.  The function must apply
the needed byte-swapping on the data in each row based on the
description of the table data in the table header.

\begin{prototype}
psArray *psFitsReadTable (psFits *fits);
\end{prototype}
This function reads the entire data block from a table into the a
\code{psArray}, with one element of the array per row.  Each row is
stored as a \code{psMetadata} collection as described above for
\code{psFitsReadTableRow}. 

\begin{prototype}
bool psFitsWriteTable(psFits* fits, const psMetadata *header, const psArray* table); 
\end{prototype}
Accepts a \code{psArray} of \code{psMetadata} and writes it to the
current HDU.  If the current HDU is not a table type, this will fail
and return FALSE.

\begin{prototype}
bool psFitsUpdateTable(psFits* fits, const psMetadata* data, int row); 
\end{prototype}
Writes the \code{psMetadata} data to a FITS table at the specified row
in the current HDU.  If the current HDU is not a table type, this will
fail and return FALSE.  

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

\pagebreak
\section{Data manipulation}

We require a variety of basic data manipulation functions which will
act upon data (in particular, arrays/vectors).  We require the
following capabilities:
\begin{itemize}
\item Vector and image arithmetic;
\item Sorting;
\item Statistics;
\item Matrix operations and linear algebra;
\item (Fast) Fourier Transforms;
\item General mathematical functions; and
\item Minimization and fitting routines.
\end{itemize}

\subsection{Sorting}

We require the ability to sort a vector.  The following function
returns the vector, sorted from the smallest (i.e.\ most negative)
value in the first element, and the largest (i.e.\ most positive)
value in the last element.  The input vector, \code{in}, may be sorted
in-place if it is also specified as the \code{out} vector. This
function is specified for input types \code{psS8, psU16, psF32,
psF64}.  The input and output vectors must have the same type.

\begin{prototype}
psVector *psVectorSort(psVector *outVector, const psVector *inVector);
\end{prototype}

We also require the ability to sort one vector based on another.  For
example, we may want to sort both \code{x} and \code{y} by the value
in \code{x}.  In order to facilitate this, we will have a sort
function return a vector containing the indices for the unsorted list
in the order appropriate for the sorted vector, sorted from the
smallest (i.e.\ most negative) value in the first element, and the
largest (i.e.\ most positive) value in the last element.  The output
vector must be of type \code{psU32}.  This function is specified for
input types \code{psS8, psU16, psF32, psF64}.

\begin{prototype}
psVector *psVectorSortIndex(psVector *outVector, const psVector *inVector);
\end{prototype}

The sorted vectors may be accessed in the following manner:
\begin{verbatim}
indexVector = psVectorSortIndex(NULL, x);
for (int i = 0; i < indexVector.n; i++) {
    doMyFunc(x[indexVector.arr.arr_U32[i]], y[indexVector[i].arr.arr_U32]);
}
\end{verbatim}


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

\subsection{Masks}

Several functions --- especially the statistical and image functions
--- use a mask vector and a mask value to specify values in an input
list that should be excluded from the calculations.  We define a mask
type, which we will initially set to U8.  In the event that our masks
exceed eight bits, we can then change the mask definition without
major changes throughout the code.

\begin{datatype}
typedef psU8 psMaskType;
\end{datatype}

\subsection{Statistical Functions}

\subsubsection{Statistical measures}

We require a very general statistics function, which, given a vector
of floating-point values, will be able to calculate the following
population statistics:
\begin{itemize}
\item Minimum value;
\item Maximum value;
\item Sample mean;
\item Sample median;
\item Sample standard deviation;
\item Sample upper and lower quartiles;
\item Clipped mean and number of values used to calculate;
\item Clipped standard deviation; 
\item Robust median and number of values used to calculate;
\item Robust standard deviation;
\item Robust upper and lower quartiles;
\item Fitted mean and number of values used to calculate;
\item Fitted standard deviation;
\end{itemize}

For definitions of each of these, see the accompanying Algorithms
Definition Document (ADD), but in general, ``sample'' refers to the
entire vector, ``clipped'' refers to clipping the distribution,
``robust'' refers to determining the quantities from the cumulative
histogram, and ``fitted'' refers to fitting the data histogram with a
Gaussian model.  Each of these must be available from a single
function:

\begin{prototype}
psStats *psVectorStats(psStats *stats,
                       const psVector *in, 
                       const psVector *errors,
                       const psVector *mask,
                       psMaskType maskVal);
\end{prototype}
%
This function takes the input data in \code{in} (with optional
\code{errors} in these values; and with optional masking in
\code{mask}, so that the user may explicitly reject specific entries)
and a \code{psStats} structure, which will be altered and returned.
The \code{psStats} structure includes several fields which are used
for both input and output: \code{nSubsample} (default value of
100,000) specifies the maximum number of data points to be used for
the statistics calculation.  If more points than this are available,
then the input vector must be randomly sampled to provide
\code{nSubsample} measurements.  \code{min} and \code{max} may be used
to specify a value range for which the statistics are calculated.
\code{binsize} specifies a choice for the robust statistics histogram
bin size.  If these are to be used, the user must set the
corresponding \code{options} bits \code{PS_STAT_USE_RANGE} or
\code{PS_STAT_USE_BINSIZE}.  \code{clipSigma} specifies the number of
standard deviations for which data should be clipped.  \code{clipIter}
specifies the number of iterations which should be used for clipping.
The defaults for these two numbers is both 3.  Since the sample
statistics scale like $N\log N$, for large numbers of input data
points, it is faster to use the robust statistics.  Default input
field values must be set by the \code{psStats} constructor.  The input
vector may be of type \code{psS8}, \code{psU16}, \code{psF32},
\code{psF64}.  The \code{errors} must be of the same type as the
\code{in} vector.  If \code{errors} is not \code{NULL}, the
calculation of certain statistics are modified: The sample mean is
calculated using the formula for the weighted mean; the standard
deviation is modified as specified in the ADD; the clipping used for
clipped statistics are modified according to the ADD; the robust
median and quartiles are modified as specified in the ADD.  The mask
must be of type \code{psMaskType}.

The \code{psStats} structure is defined with entries for each of the
desired statistical quantities:

\begin{datatype}
typedef struct {
    double sampleMean;                  ///< formal mean of sample
    double sampleMedian;                ///< formal median of sample
    double sampleStdev;                 ///< standard deviation of sample
    double sampleUQ;                    ///< upper quartile of sample
    double sampleLQ;                    ///< lower quartile of sample
    double robustMedian;                ///< robust median of data
    double robustStdev;                 ///< robust standard deviation of data
    double robustUQ;                    ///< robust upper quartile
    double robustLQ;                    ///< robust lower quartile
    int    robustN50;                   ///< Number of points UQ-LQ
    double fittedMean;                  ///< robust mean of data
    double fittedStdev;                 ///< robust standard deviation of data
    int    fittedNfit;                  ///< Number of points in Gauss. fit
    double clippedMean;                 ///< Nsigma clipped mean
    double clippedStdev;                ///< standard deviation after clipping
    int    clippedNvalues;              ///< number of data points used for clipped mean
    double clipSigma;                   ///< Nsigma used for clipping; user input
    int    clipIter;                    ///< Number of clipping iterations; user input
    double min;                         ///< minimum data value in data; input range
    double max;                         ///< maximum data value in data; input range
    double binsize;                     ///< binsize for robust fit (input/output)
    int    nSubsample;                  ///< maximum number of measuremenst (input)
    psStatsOptions options;             ///< bitmask of calculated values
} psStats;
\end{datatype}
where \code{psStatsOptions} is defined with entries to turn on the
calculation of each of the statistics:

\begin{datatype}
typedef enum {
    PS_STAT_SAMPLE_MEAN           = 0x000001,
    PS_STAT_SAMPLE_MEDIAN         = 0x000002,
    PS_STAT_SAMPLE_STDEV          = 0x000004,
    PS_STAT_SAMPLE_QUARTILE       = 0x000008, 
    PS_STAT_ROBUST_MEDIAN         = 0x000010,
    PS_STAT_ROBUST_STDEV          = 0x000020,
    PS_STAT_ROBUST_QUARTILE       = 0x000040, 
    PS_STAT_FITTED_MEAN           = 0x000080,
    PS_STAT_FITTED_STDEV          = 0x000100,
    PS_STAT_CLIPPED_MEAN          = 0x000200,
    PS_STAT_CLIPPED_STDEV         = 0x000400,
    PS_STAT_MAX                   = 0x000800,     
    PS_STAT_MIN                   = 0x001000,
    PS_STAT_USE_RANGE             = 0x002000,
    PS_STAT_USE_BINSIZE           = 0x004000
} psStatsOptions;                         
\end{datatype}

A constructor is also required:
%
\begin{prototype}
psStats *psStatsAlloc(psStatsOptions options);
\end{prototype}

\begin{prototype}
long psVectorCountPixelMask (psVector *mask, psMaskValue value);
\end{prototype}
This function returns the number of pixels in the \code{vector} which
satisfy any of the mask bits.  An error (eg, invalid vector) results
in a return value of -1.  The vector must be of the same type as
psMaskValue.

\subsubsection{Histograms}
\label{sec:histograms}

We also require to be able to generate histograms, given a list of
upper and lower bounds for each of the bins.  We define the following
data structure to represent a histogram:
\begin{datatype}
typedef struct {
    const psVector *bounds;             ///< Bounds for the bins
    psVector *nums;                     ///< Number in each of the bins
    int minNum;                         ///< Number below minimum
    int maxNum;                         ///< Number above maximum
    bool uniform;                       ///< Is it a uniform distribution?
} psHistogram;
\end{datatype}
In this structure, the vector \code{bounds} specifies the boundaries
of the histogram bins, and must of type \code{psF32}, while
\code{nums} specifies the number of entries in the bin, and must be of
type \code{psF32} in order to accomodate errors.  The value of
\code{bounds.n} must therefore be 1 greater than than \code{nums.n}.
The two values \code{minNum} and \code{maxNum} are the number of data
values which fell below the lower limit bound or above the upper limit
bound, respectively.

The constructors and destructor follow.  We specify two constructors,
so that the bounds of the bins may either be specified explicitly, or
implicitly through simply specifying an upper and lower limit along
with the size of the bins.

\begin{prototype}
psHistogram *psHistogramAlloc(float lower, float upper, int n);
\end{prototype}
where \code{lower} specifies the lower bound of the histogram range,
\code{upper} specified the upper bound of the histogram range, and
\code{n} is the number of bins desired across the range.  This
constructor sets the value of \code{uniform} to be true.

A histogram with a more flexible bin set may be constructed with the
following constructor:
\begin{prototype}
psHistogram *psHistogramAllocGeneric(const psVector *bounds);
\end{prototype}
where the \code{psVector *bounds} (of type \code{psF32}) is defined by
the user to specify the boundaries of the histogram bins. This
constructor sets the value of \code{uniform} to be false.

The following function populates the histogram bins from the specified
vector (\code{values}), and optionally the \code{errors} in the input
values.  It alters and returns the histogram \code{out} structure.
\begin{prototype}
psHistogram *psVectorHistogram(psHistogram *out, 
                               const psVector *values,
                               const psVector *errors,
                               const psVector *mask,
                               psMaskType maskVal);
\end{prototype}
The \code{values} vector may be of types \code{psU8, psU16, psF32,
psF64}, with the \code{errors} vector of the corresponding type.

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

\subsection{Analytical functions}

\subsubsection{Polynomials}

PSLib provides APIs to represent and interact with polynomials in up
to four dimensions, with both floating-point and double-precision
numbers.  In Pan-STARRS processing, the astrometry requirements push
the need for at least four dimensions ($x$,$y$, color and magnitude)
and double-precision (for milli-arcsec precision) versions.  We must
also be able to calculate the errors in the fit coefficients, as well
as be able to turn on and off each coefficient.

In addition to general polynomials ($\sum_{i=0}^n a_i x^i$), we also
use Chebyshev polynomials ($\sum_{i=0}^n a_i T_i(x)$), which have
properties which are useful in the modeling of data over a defined
domain.  Note that Chebyshev polynomials should only have inputs in
the range $-1 \le x \le +1$ (because they are bounded over this
range), but we will not enforce this.

This leads us to define the following polynomial types:

\begin{datatype}
/** One-dimensional polynomial */
typedef struct {
    psPolynomialType type;              ///< Polynomial type
    unsigned int n;                     ///< Number of terms
    psF64 *coeff;                       ///< Coefficients
    psF64 *coeffErr;                    ///< Error in coefficients
    char *mask;                         ///< Coefficient mask
} psPolynomial1D;
\end{datatype}

\begin{datatype}
/** Two-dimensional polynomial */
typedef struct {
    psPolynomialType type;              ///< Polynomial type
    unsigned int nX                     ///< Number of terms in x
    unsigned int nY;                    ///< Number of terms in y
    psF64 **coeff;                      ///< Coefficients
    psF64 **coeffErr;                   ///< Error in coefficients
    char **mask;                        ///< Coefficients mask
} psPolynomial2D;
\end{datatype}

etc., up to four dimensions.  \code{psPolynomialType} is an
enumerated type specifying the type of the polynomial: ordinary
or Chebyshev:

\begin{datatype}
typedef enum {
    PS_POLYNOMIAL_ORD,                  ///< Ordinary polynomial
    PS_POLYNOMIAL_CHEB                  ///< Chebyshev polynomial
} psPolynomialType;
\end{datatype}

The constructors are:
\begin{prototype}
psPolynomial1D *psPolynomial1DAlloc(psPolynomialType type, 
                                    unsigned int nX);
psPolynomial2D *psPolynomial2DAlloc(psPolynomialType type, 
                                    unsigned int nX, unsigned int nY);
psPolynomial3D *psPolynomial3DAlloc(psPolynomialType type, 
                                    unsigned int nX, unsigned int nY, unsigned int nZ);
psPolynomial4D *psPolynomial4DAlloc(psPolynomialType type, 
                                    unsigned int nX, unsigned int nY, unsigned int nZ, unsigned int nT);
\end{prototype}
where \code{nX}, \code{nY}, etc specify the polynomial order in the
given dimension.  The coefficients, errors and masks are set initially
to zero.

To evaluate the polynomials at specific coordinates, we define:
\begin{prototype}
psF64 psPolynomial1DEval(const psPolynomial1D *poly, 
                         psF64 x);
psF64 psPolynomial2DEval(const psPolynomial2D *poly, 
                         psF64 x, 
                         psF64 y);                
psF64 psPolynomial3DEval(const psPolynomial3D *poly, 
                         psF64 x, 
                         psF64 y, 
                         psF64 z);                 
psF64 psPolynomial4DEval(const psPolynomial4D *poly, 
                         psF64 x, 
                         psF64 y, 
                         psF64 z, 
                         psF64 t);
\end{prototype}

In the event that several evaluations are required, we also define:
\begin{prototype}
psVector *psPolynomial1DEvalVector(const psPolynomial1D *poly, 
                                   const psVector *x);
psVector *psPolynomial2DEvalVector(const psPolynomial2D *poly, 
                                   const psVector *x, 
                                   const psVector *y);
psVector *psPolynomial3DEvalVector(const psPolynomial3D *poly, 
                                   const psVector *x, 
                                   const psVector *y, 
                                   const psVector *z);
psVector *psPolynomial4DEvalVector(const psPolynomial4D *poly, 
                                   const psVector *x, 
                                   const psVector *y, 
                                   const psVector *z, 
                                   const psVector *t);
\end{prototype}
The function shall accept input vectors of any type, converting to
\code{psF64} as needed.  In the event that the \code{x} and \code{y}
vectors are of differing sizes, the function shall generate a warning,
truncate the longer vector to the size of the shorter, and continuing.
The precision of the output \code{psVector} shall be \code{psF64}.  In
evaluation, those coefficients that have the corresponding \code{mask}
element non-zero shall not be evaluated.

\subsubsection{Splines}

A spline is a popular choice for fitting 1D data, such as overscans,
but we neglected to define them for PSLib.  We now define
one-dimensional cubic splines, \code{psSpline1D}, which shall be
incorporated into PSLib:

\begin{datatype}
typedef struct {
    unsigned int n;                     ///< Number of spline pieces
    psPolynomial1D **spline;            ///< Array of n pointers to splines
    psVector *knots;                    ///< The boundaries between each spline piece.  Size is n+1.
} psSpline1D;
\end{datatype}

The \code{psSpline1D} structure consists of an array of \code{n}
polynomials, which are the spline pieces.  Note that this means that
the spline pieces may, in general, be of any order.  \textbf{For the
present, we shall restrict the order of the polynomials to either 1
(linear) or 3 (cubic).}  All the spline pieces shall have the same
order polynomial (the type of polynomial is left to the
implementation).  The \code{knots} member specifies the boundaries
between each spline piece (including the two ends).  The \code{knots}
vector may be of type U32 or F32.

Of course, we require the appropriate constructors and destructor:
\begin{prototype}
psSpline1D *psSpline1DAlloc(unsigned int n, unsigned int order, float min, float max);
psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, unsigned int order);
\end{prototype}

\code{psSpline1DAlloc} shall allocate and return a \code{psSpline1D},
setting the \code{knots} on the basis of the input number of spline
pieces, \code{n}, and the minimum (\code{min}) and maximum
(\code{max}) data values.  The spline pieces shall be of the specified
\code{order}.

\code{psSpline1DAllocGeneric} shall allocate and return a
\code{psSpline1D}, using the \code{bounds} to define the number of
spline pieces and the \code{knots}.  The spline pieces shall be of
the specified \code{order}.

Also, as for the polynomials, we require evaluators.  Given a
\code{spline} and ordinate at which to evaluate, \code{x},
\code{psSpline1DEval} shall evaluate and return the value of the
spline at the ordinate.  If the ordinate is outside the bounds, then
the function shall generate a warning, and extrapolate the spline to
the ordinate and return the value.  Similarly,
\code{psSpline1DEvalVector} shall return a vector of evaluated values
for an input vector of ordinates.

\begin{prototype}
float psSpline1DEval(const psSpline1D *spline, float x);
psVector *psSpline1DEvalVector(const psSpline1D *spline, const psVector *x);
\end{prototype}

\subsubsection{Gaussians}

Gaussians are used extensively in any data-analysis system, in
particular to represent a real population distribution.  We require a
function to evaluate a Gaussian for a given coordinate.

The Gaussian evaluation is provide by:
\begin{prototype}
float psGaussian(float x, float mean, float sigma, bool normal);
\end{prototype}
which evaluates a Gaussian with the given \code{mean} and \code{sigma}
at the given coordinate \code{x}.  If \code{normal} is true, the
Gaussian is normalized (total integral = 1), otherwise, the Gaussian
is non-normalized (central peak value = 1).  The evaluated Gaussian
is:

\[ \frac{1}{\sqrt{2\pi\sigma^2}} exp^{-\frac{(x-mean)^2}{2\sigma^2}} \]

In the case of the non-normalized Gaussian, the leading coefficient is
dropped.

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

\subsection{Minimization and fitting routines}

We require a general minimization routine, a routine that will
specifically minimize $\chi^2$ given a list of data with associated
errors, and functions that will analytically determine the best
polynomial and spline fits by minimizing $\chi^2$.

We specify two minimization engines, Levenberg-Marquardt and Powell,
since the former is relatively robust, but requires derivatives to be
known, while the latter does not require derivatives to be known, but
since it needs to calculate the derivatives on the fly, is more
computationally intensive.  

We define the structure \code{psMinimization} to carry in user-defined
limits on the minimization process, and to carry out information about
the minimization analysis.  The maximum number of iterations is
specified by \code{maxIter}, while the maximum tolerance for
convergence is \code{tol}.  Four parameter vectors describe the
behavior of individual parameters.  The \code{paramMask} vector
defines the free (0) or frozen (not 0) parameters.  The
\code{paramMin} defines the minimum allowed value for each parameter,
while \code{paramMax} defines the maximum allowed value for each
parameter.  During the analysis, parameters which would trend beyond
these limits saturate to the limit.  The \code{paramDelta} specifies
the maximum allowed absolute value of the swing of a given parameter.
If the delta for a specific parameter would be larger than this, the
swing is saturated to the limit (with the correct sign).  Any of these
parameter vectors may be set to \code{NULL}, in which case the concept
is ignored in the analysis.  The output information carried by the
structure consists of the value of the function at the minimum
(\code{value}), the number of iterations performed (\code{iter}) and
last change in tolerance before returning (\code{lastDelta}).

\begin{datatype}
typedef struct {
    const int maxIter;                  ///< Maximum number of iterations
    const float tol;                    ///< Tolerance to reach
    float value;                        ///< Value after minimization
    int iter;                           ///< Actual number of iterations performed
    float lastDelta;                    ///< Last change before quitting
    psVector *paramMask;                ///< valid / invalid parameters
    psVector *paramMax;                 ///< max allowed parameters
    psVector *paramMin;                 ///< min allowed parameters
    psVector *paramDelta;               ///< max allowed param swing
} psMinimization;
\end{datatype}

The corresponding allocator is:
\begin{prototype}
psMinimization *psMinimizationAlloc(int maxIter, float tol);
\end{prototype}
and the parameter vectors are initially set to \code{NULL}.  

\subsubsection{Levenberg-Marquardt}

Consider a function of a collection of parameters, \code{params},
which is evaluated at a position, \code{x}, which returns a single
floating point value which is the value of the function given the
parameters and coordinate vectors, along with the derivatives of the
function with respect to each of the parameters, \code{deriv}:
\begin{datatype}
typedef float (*psMinimizeLMChi2Func)(psVector *deriv, const psVector *params, const psVector *x);
\end{datatype}

Then \code{psMinimizeLMChi2} shall fit the specified function,
\code{func}, to a set of measurements, \code{x,y,yErr}, using the
Levenberg-Marquardt method:

\begin{prototype}
bool psMinimizeLMChi2(psMinimization *min, 
                      psImage *covar, 
                      psVector *params,
                      const psArray  *x, 
                      const psVector *y,
                      const psVector *yErr, 
                      psMinimizeLMChi2Func func);
\end{prototype}

The function shall return \code{false} in the event that there was an
error (bad input, too many iterations), and also call \code{psError}.
It is not an error for the minimization to reach one of the parameter
limits.  

The minimization specification, \code{min}, shall be modified with the
\code{value}, \code{iter} and \code{lastDelta} members updated with
the values appropriate from the minimization.

On calling the minimizer, the \code{params} vector shall consist of a
``best guess'' for the parameters that minimize the model function,
\code{func}.  On successful completion, the input parameters,
\code{params}, shall be updated with the values that minimize the
input function.  The function shall also update the covariance matrix,
\code{covar}, in the event that it is non-\code{NULL}.  The function
shall generate in error in the event that \code{covar} is
non-\code{NULL} and is not a square matrix with size matching that of
\code{params}.

Parameters that have a corresponding \code{paramMask} entry that is
non-zero are to be held fixed by the minimizer.  It shall be an error
for \code{paramMask} not to be of the same dimension as \code{params}.

The measurement ordinates, \code{x}, shall consist of multiple
vectors, each of which may be passed to the model \code{func}.  If the
measurement coordinates, \code{y}, and errors, \code{yErr}, are not of
the same length as the ordinates array, \code{x}, then the function
shall generate a warning, and truncate the longest of the
array/vectors to match the length of the shortest.  The vectors
contained within the \code{x} array, and the \code{y} and \code{yErr}
vectors must be of type \code{psF32}.  The \code{yErr} vector may be
\code{NULL}, in which case the errors shall be assumed to be
identical.

\code{paramMask} must be of type \code{psMaskType}, while \code{params}
must be of type \code{psF32}.  The \code{func} function must be valid
only for types \code{psF32}, \code{psF64}.

\begin{prototype}
bool psMinimizeGaussNewtonDelta(psVector *delta, 
                                const psVector *params,
                                const psVector *paramMask, 
                                const psArray  *x, 
                                const psVector *y,
                                const psVector *yErr, 
                                psMinimizeLMChi2Func func);
\end{prototype}

The function \code{psMinimizeGaussNewtonDelta} can be used to evaluate
the goodness of a fit.  This function returns the distances of the
current parameter set from the minimization parameters expected from
Gauss-Newton minimization.  The arguments have the same meaning as for
\code{psMinimizeLMChi2}.  The returned vector, \code{delta}, consists
of the distances.  This vector must be pre-allocated to the
dimenstions of \code{params}.

%% \subsubsubsection{Pre-defined Functions for LM}
%% 
%% We define some commonly used functions for use with the LM
%% minimization, used for the purpose of performing $\chi^2$ fitting:
%% 
%% \begin{prototype}
%% psMinimizeLMChi2Func psMinimizeLMChi2Gauss1D;
%% psMinimizeLMChi2Func psMinimizeLMChi2Gauss2D;
%% \end{prototype}
%% 
%% \code{psMinimizeChi2LMGauss1D} shall take as \code{params}, the
%% normalization, center, and standard deviation of a Gaussian to be fit,
%% and as \code{x}, a vector containing a single value.  It shall return
%% the value of the Gaussian at the value, and the derivatives
%% (\code{deriv}) with respect to each of the parameters.
%% 
%% \code{psMinimizeChi2LMGauss2D} shall take, as \code{params}, the
%% normalization, center (two values), standard deviation (two values)
%% and position angle of a 2-dimensional Gaussian, and as \code{x}, a
%% vector containing a position, $(x,y)$.  It shall return the value of
%% the 2-dimensional Gaussian at the specified point, along with the
%% derivatives with respect to each of the parameters.

\subsubsection{Powell}

As for the LM minimizer, consider a function of a collection of
parameters, \code{params}, and (possibly several) coordinate vectors
(which we represent as an array of vectors), \code{coord}, and returns
a single floating point value which is the value of the function given
the parameters and coordinate vectors, but now the derivatives are not
known:
\begin{datatype}
typedef float (*psMinimizePowellFunc)(const psVector *params, const psArray *coords);
\end{datatype}

Then \code{psMinimizePowell} shall minimize the specified function,
\code{func}, using the Powell method:

\begin{prototype}
bool psMinimizePowell(psMinimization *min, 
                      psVector *params, 
                      const psVector *paramMask,
                      const psArray *coords, 
                      psMinimizePowellFunc func);
\end{prototype}

The inputs and general behavior of this function is the same as for
\code{psMinimizeLM}, except for the absence of the covariance matrix,
\code{covar} \note{why does this not have a covar matrix?}

\subsubsubsection{Minimizing $\chi^2$ with Powell}

We require a front-end to the Powell minimizer to allow fitting
general functions to data.

\begin{datatype}
typedef psVector* (*psMinimizeChi2PowellFunc)(const psVector *params, const psArray *coords);
\end{datatype}

\begin{prototype}
bool psMinimizeChi2Powell(psMinimization *min, 
                          psVector *params, 
			  const psVector *paramMask,
                          const psArray *coords, 
			  const psVector *value, 
			  const psVector *error,
                          psMinimizeChi2PowellFunc model);
\end{prototype}

The inputs and general behavior of \code{psMinimizeChi2Powell} is
similar as for \code{psMinimizePowell}, with the exception that
instead of being provided the function to be minimized, it is provided
a model to fit and must calculate the $\chi^2$ to be minimized itself.

\note{why does this not have a covar matrix?}

\note{unify the param names with psMinimizeLMChi2Func?}

For this purpose, \code{value} shall contain measured values at the
coordinates, and \code{error} may either be non-\code{NULL}, in which
case it contains the errors in the measured values; otherwise the
errors shall assumed to be unity for the purpose of fitting where the
errors are all identical (and possibly unknown).

Furthermore, the \code{model} function provided by the user shall
return a vector of values (instead of a single value, as was the case
for \code{psMinimizePowell}), corresponding to the model predictions
to be compared against the measured values.  If the lengths of the
\code{value}, \code{error} (if provided) vectors and the vector
returned by the \code{model} function differ, then
\code{psMinimizeChi2Powell} shall generate a warning, truncate the
vectors to the length of the shortest and proceed (only one error
should be generated per call).

\subsubsection{Analytical fits}

\begin{prototype}
psPolynomial1D *psVectorFitPolynomial1D(psPolynomial1D *poly, 
                                        const psVector *mask, 
                                        psMaskType      maskValue,
                                        const psVector *f,
                                        const psVector *fErr,
                                        const psVector *x);
psPolynomial2D *psVectorFitPolynomial2D(psPolynomial2D *poly, 
                                        const psVector *mask, 
                                        psMaskType      maskValue,
                                        const psVector *f,
                                        const psVector *fErr,
                                        const psVector *x, 
                                        const psVector *y);
psPolynomial3D *psVectorFitPolynomial3D(psPolynomial3D *poly, 
                                        const psVector *mask, 
                                        psMaskType      maskValue,
                                        const psVector *f,
                                        const psVector *fErr,
                                        const psVector *x, 
                                        const psVector *y, 
                                        const psVector *z);
psPolynomial4D *psVectorFitPolynomial4D(psPolynomial4D *poly, 
                                        const psVector *mask, 
                                        psMaskType      maskValue,
                                        const psVector *f,
                                        const psVector *fErr,
                                        const psVector *x, 
                                        const psVector *y, 
                                        const psVector *z, 
                                        const psVector *t);
\end{prototype}
These functions return the polynomial that best fits the input data.
The provided polynomial, \code{poly}, specifies the fit order, and
will be returnd with the best-fit coefficients.  The dependent
variable and its error (\code{f, fErr}) are provided along with the
appropriate number of independent variables (\code{x}, \code{y}, etc).
In the special case of 1D fitting, the independent variable list,
\code{x} may be \code{NULL}, in which case the vector index is used.
The dependent variable error, \code{yErr} may be null, in which case
the solution is determined in the assumption that all data errors are
equal.  This function must be valid only for input data types of
\code{psF32}, \code{psF64}. The \code{mask} and \code{maskValue}
entries may be used to specify an input data set to ignore.  All of
the input vectors must be the same length.

\begin{prototype}
psPolynomial1D *psVectorClipFitPolynomial1D(psPolynomial1D *poly, 
                                        psStats            *stats,
                                        const psVector     *mask, 
                                        psMaskType          maskValue,
                                        const psVector     *f,
                                        const psVector     *fErr,
                                        const psVector     *x);
psPolynomial2D *psVectorClipFitPolynomial2D(psPolynomial2D *poly, 
                                        psStats            *stats,
                                        const psVector     *mask, 
                                        psMaskType          maskValue,
                                        const psVector     *f,
                                        const psVector     *fErr,
                                        const psVector     *x, 
                                        const psVector     *y);
psPolynomial3D *psVectorClipFitPolynomial3D(psPolynomial3D *poly, 
                                        psStats            *stats,
                                        const psVector     *mask, 
                                        psMaskType          maskValue,
                                        const psVector     *f,
                                        const psVector     *fErr,
                                        const psVector     *x, 
                                        const psVector     *y, 
                                        const psVector     *z);
psPolynomial4D *psVectorClipFitPolynomial4D(psPolynomial4D *poly, 
                                        psStats            *stats,
                                        const psVector     *mask, 
                                        psMaskType          maskValue,
                                        const psVector     *f,
                                        const psVector     *fErr,
                                        const psVector     *x, 
                                        const psVector     *y, 
                                        const psVector     *z, 
                                        const psVector     *t);
\end{prototype}
These functions replicate the capability of the vector fitting
functions, but also include an iterative clipping stage.  The clipping
parameters are defined by the \code{stats} entry, to which are also
returned the clipped statistics for the final residuals.  Thus, if N
clipping iterations are requested, the function performs the fit,
constructs the residuals, measures the residual scatter, and masks the
outlier vector elements.  This cycle is repeated N times, though on
the last iteration, no additional masking is performed.  The provided
mask must be respected, and any additionally masked elements are also
masked with the same mask vector, which must be provided.  The
elements masked by this routine are given the mask value of 1 and may
thus be distinguished from input masked elements if they use a
different mask value.

\begin{prototype}
 bool psVectorFitSpline1D(psSpline1D *spline, const psVector *x, const psVector *y,
                          const psVector * yErr);
\end{prototype}
\code{psVectorFitSpline1D} shall fit the code{spline} that best fits
the given combination of ordinates (\code{x}) and coordinates
(\code{y}) with the known errors (\code{yErr}, which may be
\code{NULL}).  As is the case for \code{psVectorFitPolynomial1D}, if
\code{x} is \code{NULL}, then the index of \code{y} shall be used as
the ordinate.  This function must be valid only for types
\code{psF32}, \code{psF64}.

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

\subsection{Image Operations}

We require a variety of functions to manipulate these image
structures, including creation, destruction, input, output, and
various manipulations of the pixels.  The required functions are
listed below, and fall into several categories.

\subsubsection{Image Structure Manipulation}

\begin{prototype}
psImage *psImageSubset(psImage *image, psRegion region);
\end{prototype}
Define a subimage of the specified area of the given image.  This
function must raise an error if the requested subset area lies outside
of the parent image and return \code{NULL}.  The argument \code{image}
is the parent image, \code{region.x0, region.y0} specify the starting
pixel of the subraster, and \code{region.x1,region.y1} specify the
extent of the desired subraster.  Note that the row and column of this
``upper right-hand corner'' are \textit{NOT} included in the region.
In the event that \code{x1} or \code{y1} are negative, they shall be
interpreted as being relative to the size of the parent image in that
dimension.  The entire subraster must be contained within the raster
of the parent image.  Note that the \code{refCounter} for the parent
should be incremented.  This function must be defined for the
following types: \code{psU8}, \code{psU16}, \code{psS8}, \code{psS16},
\code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.

\begin{prototype}
psImage *psImageCopy(psImage *output, const psImage *input, psElemType type);
\end{prototype}
Create a copy of the specified image, converting the type in the
process.  If the output target pointer is not \code{NULL}, place the
result in the specified structure.  The output image data must be
allocated as a single, contiguous block of memory.  The output image
may not be the input image.  The \code{col0,row0} of the \code{input}
image shall be preserved.  This function must be defined for the
following types: \code{psU8}, \code{psU16}, \code{psS8}, \code{psS16},
\code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.

\begin{prototype}
psImage *psImageTrim(psImage *image, psRegion region);
\end{prototype}
Trim the specified \code{image} in-place, which involves shuffling the
pixels around in memory.  The region to be kept is defined by the
lower-left corner, \code{region.x0,region.y0}, and the upper-right
corner, \code{region.x1,region.y1}.  Note that the row and column of
the ``upper right-hand corner'' are \textit{NOT} included in the
region.  In the event that \code{region.x1} or \code{region.y1} are
negative, they shall be interpreted as being relative to the size of
the parent image in that dimension.

The function shall generate an error if the specified region is
outside the bounds of the input image.  Any children of the input
image shall be freed by \code{psImageTrim} before the trim takes
place.

\subsubsection{Image Pixel Extractions}

\begin{datatype}
typedef enum {
    PS_CUT_X_POS,                     ///< Cut in positive x direction
    PS_CUT_X_NEG,                     ///< Cut in negative x direction
    PS_CUT_Y_POS,                     ///< Cut in positive y direction
    PS_CUT_Y_NEG                      ///< Cut in negative y direction
} psImageCutDirection;
\end{datatype}

\begin{prototype}
psVector *psImageRow(psVector *out, 
                     const psImage *input,
		     psU32 row);
psVector *psImageCol(psVector *out, 
                     const psImage *input,
		     psU32 column);
\end{prototype}

These two functions extract a single complete \code{row} or
\code{column} from the \code{image} and return it to the provided
\code{vector}, allocating it if it is \code{NULL}.  These are provided
as fast, simple implementations of data extraction.  For more
sophisticated operations, the following two functions should be used. 

\begin{prototype}
psVector *psImageSlice(psVector *out, 
                       psPixels *coords,
                       const psImage *input,
                       const psImage *mask, 
                       psMaskType maskVal, 
                       psRegion region,
                       psImageCutDirection direction, 
                       const psStats *stats);
\end{prototype}
Extract pixels from rectilinear region to a vector (array of floats).
The output vector contains either \code{region.x1-region.x0} or
\code{region.y1-region.y0} elements, based on the value of the
direction: e.g., if \code{direction} is \code{PS_CUT_X_POS}, there are
\code{region.x1-region.x0} elements.  The region to be ``sliced'' is
defined by the lower-left corner, \code{region.x0,region.y0}, and the
upper-right corner, \code{region.x1,region.y1}.  Note that the row and
column of the ``upper right-hand corner'' are \textit{NOT} included in
the region.  In the event that \code{region.x1} or \code{region.y1}
are negative, they shall be interpreted as being relative to the size
of the parent image in that dimension.

The input region is collapsed in the direction perpendicular to that
specified by \code{direction}, and each element of the output vectors
is derived from the statistics of the pixels at that direction
coordinate.  The statistic used to derive the output vector value is
specified by \code{stats}.  If \code{mask} is non-\code{NULL}, pixels
for which the corresponding \code{mask} pixel matches \code{maskVal}
are excluded from operations.  If \code{coords} is not \code{NULL},
the calculated coordinates along the slice are returned in this
array of pixels.  Only one of the statistics choices may be specified,
otherwise the function must return an error.  This function must be
defined for the following types: \code{psS8}, \code{psU16},
\code{psF32}, \code{psF64}.

\begin{prototype}
psVector *psImageCut(psVector *out, 
                     psVector *cutCols, 
                     psVector *cutRows, 
                     const psImage *input, 
                     const psImage *mask, 
                     psMaskType maskVal,
                     psRegion region,
                     unsigned int nSamples, 
                     psImageInterpolateMode mode);
\end{prototype}
Extract pixels along a line segment, \code{(region.x0,region.y0)} to
\code{(region.x1,region.y1)}, on the image to a vector of the same
data type.  The line segment is sampled \code{nSamples} times, hence
the output vector contains \code{nSamples} elements.  The
interpolation method used to derive the output vector value at each
sample along the line segment is specified by \code{mode}.  If the
\code{mask} is non-\code{NULL}, then pixels for which the
corresponding mask value is \code{maskVal} are not included in the
interpolation.  This function must be defined for the following types:
\code{psS8}, \code{psU16}, \code{psF32}, \code{psF64}.

\begin{prototype}
psVector *psImageRadialCut(psVector *out, 
                           const psImage *input, 
                           const psImage *mask, 
                           psMaskType maskVal,
                           float x, 
                           float y,
                           const psVector *radii, 
                           const psStats *stats);
\end{prototype}
Extract radial region data to a vector.  A vector is constructed where
each vector elements is derived from the statistics of the pixels
which land within one of a sequence of radii.  The radii are centered
on the image pixel coordinate \code{x,y}, and are defined by the
sequence of values in the vector \code{radii}.  The specific algorithm
which must be used is described in the PSLib ADD (PSDC-430-006).  The
statistic used to derive the output vector value is specified by
\code{stats}.  Only one of the statistics choices may be specified,
otherwise the function must return an error.  If \code{mask} is
non-\code{NULL}, pixels for which the corresponding \code{mask} pixel
matches \code{maskVal} are excluded from operations.  This function
must be defined for the following types: \code{psS8}, \code{psU16},
\code{psF32}, \code{psF64}.

\subsubsection{Image Geometry Manipulation}

Several functions which manipulate the image geometry require the
specification of the interpolation scheme to be used.  This
information is carried by the following enum:
\begin{datatype}
typedef enum {
    PS_INTERPOLATE_FLAT, 
    PS_INTERPOLATE_BILINEAR, 
    PS_INTERPOLATE_LANCZOS2,
    PS_INTERPOLATE_LANCZOS3,
    PS_INTERPOLATE_LANCZOS4,
    PS_INTERPOLATE_BILINEAR_VARIANCE,
    PS_INTERPOLATE_LANCZOS2_VARIANCE,
    PS_INTERPOLATE_LANCZOS3_VARIANCE,
    PS_INTERPOLATE_LANCZOS4_VARIANCE
} psImageInterpolateMode;
\end{datatype}

The first five are fairly straightforward.  The last four, however,
require some explanation.  When attempting to account for noise in a
CCD image, it is customary to carry an additional image containing the
variance for each pixel.  When operating on the CCD image (performing
some transformation), we want to perform the same operation on the
variance image, but the weights of the different input pixels
contributing to the output pixel must be squared in order to propagate
the noise (adding in quadrature).

\begin{prototype}
psImage *psImageRebin(psImage *out, const psImage *in, 
                      const psImage *mask, 
                      psMaskType maskVal,
                      int scale, const psStats *stats);
\end{prototype}
Rebin image to new scale.  A new image is constructed in which the
dimensions are reduced by a factor of \code{1 / scale}.  The
\code{scale}, always a positive number, is equal in each dimension and
specified the number of pixels used to define a new pixel in the
output image.  The output image is generated from all input image
pixels.  Care must be taken on the image boundary if the image
dimensions are not divisible by the scaling factor.  In those regions,
the output pixel must be constructed from the available input pixels.
Each pixel in the output image is derived from the statistics of the
corresponding set of input image pixels based on the statistics
specified by \code{stats}.  Only one of the statistics choices may be
specified, otherwise the function must return an error.  If
\code{mask} is non-\code{NULL}, pixels for which the corresponding
\code{mask} pixel matches \code{maskVal} are excluded from operations.
This function must be defined for the following types: \code{psS8},
\code{psU16}, \code{psF32} and \code{psF64}.

\begin{prototype}
psImage *psImageResample(psImage *out, const psImage *in, 
                         int scale, psImageInterpolateMode mode);
\end{prototype}
Resample image to new scale.  A new image is constructed in which the
dimensions are increased by a factor of \code{scale}.  The
\code{scale}, always a positive number, is equal in each dimension.
The output image is generated from all input image pixels.  Each pixel
in the output image is derived by interpolating between neighboring
pixels using the specified interpolation method (\code{mode}).

\begin{prototype}
psImage *psImageRotate(psImage *out, const psImage *input, float angle,
                       complex exposed, psImageInterpolateMode mode);
\end{prototype}
Rotate the input image by given angle, specified in radians.  The
output image must contain all of the pixels from the input image in
their new frame.  Pixels in the output image which do not map to input
pixels should be set to \code{exposed}, cast to the same type as the
image.  The center of rotation is always the center pixel of the
image.  The rotation is specified in the sense that a positive angle
is an anti-clockwise rotation.  This function must be defined for the
following types: \code{psU8}, \code{psU16}, \code{psS8}, \code{psS16},
\code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.

\begin{prototype}
psImage *psImageShift(psImage *out, const psImage *input, 
                      float dx, float dy, complex exposed, psImageInterpolateMode mode);
\end{prototype}
Shift image by an arbitrary number of pixels (\code{dx,dy}) in either
direction.  If the shift values are fractional, the output pixel
values should interpolate between the input pixel values.  The output
image has the same dimensions as the input image.  Pixels which fall
off the edge of the output image are lost.  Newly exposed pixels are
set to the value given by \code{exposed}, cast to the same type as the
image.  This function must be defined for the following types:
\code{psU8}, \code{psU16}, \code{psS8}, \code{psS16}, \code{psF32},
\code{psF64}, \code{psC32}, \code{psC64}.

\begin{prototype}
psImage *psImageRoll(psImage *out, const psImage *input, int dx, int dy);
\end{prototype}
Roll image by an integer number of pixels (\code{dx,dy}) in either
direction.  The output image is the same dimensions as the input
image.  Edge pixels wrap to the other side (no values are lost).  This
function must be defined for the following types: \code{psU8},
\code{psU16}, \code{psS8}, \code{psS16}, \code{psF32}, \code{psF64},
\code{psC32}, \code{psC64}.

\begin{prototype}
psImage *psImageTransform(psImage *output, 
                          psArray **blankPixels, 
                          const psImage *input,
                          const psImage *inputMask, 
                          psMaskType inputMaskVal, 
                          const psPlaneTransform *outToIn,
                          psRegion region, 
                          const psPixels *pixels, 
                          psImageInterpolateMode mode,
                          double exposedValue);
\end{prototype}
Transform the \code{input} image according the supplied
transformation.  The size of the transformed image is defined by the
\code{region} (size \code{region.x1 - region.x0} by \code{region.y1 -
region.y0}, with \code{out->x0 = region.x0} and \code{out->y0 =
region.y0}).

If the \code{inputMask} is non-\code{NULL}, those pixels in the
\code{inputMask} matching \code{inputMaskVal} are to be ignored in the
transformation.  The \code{inputMask} must be of type \code{psMaskType}, and
of the same size as the \code{input}, otherwise the function shall
generate an error and return \code{NULL}.  The transformation
\code{outToIn} specifies the coordinates in the input image of a pixel
in the output image --- note that this is the reverse of what might be
naively expected, but it is what is required in order to use
\code{psImagePixelInterpolate}.  If the \code{pixels} array is
non-\code{NULL}, it shall consist of \code{psPixelCoord}s, and only
those pixels in the output image shall be transformed; otherwise, the
entire image is generated.  The interpolation is performed using the
specified interpolation \code{mode}.  Where a pixel in the output
image does not correspond to a pixel in the input image (or all
appropriate pixels in the input image are masked), the value shall be
set to \code{exposed}, and the pixel added to the appropriate list of
pixels (\code{psPixels}) in the array of \code{blankPixels} for return
to the user.  This function must be capable of handling the following
types for the \code{input} (with corresponding types for the
\code{output}): \code{psF32}, \code{psF64}.


\subsubsection{Image Statistical Functions}

\begin{prototype}
psStats *psImageStats(psStats *stats, 
                      const psImage *in,
                      const psImage *mask,
                      psMaskType maskVal);
\end{prototype}
Determine statistics for image (or subimage).  The statistics to be
determined are specified by \code{stats}.  The \code{mask} allows
pixels to be excluded if their corresponding mask pixel value matches
the value of \code{maskVal}.  This function must be defined for the
following types: \code{psS8}, \code{psU16}, \code{psF32},
\code{psF64}.

\begin{prototype}
psHistogram *psImageHistogram(psHistogram *out, 
                              const psImage *in,
                              const psImage *mask, 
                              psMaskType maskVal);
\end{prototype}
Construct a histogram from an image (or subimage).  The histogram to
generate is specified by \code{psHistogram hist} (see
section~\ref{sec:histograms}).  The \code{mask} and \code{maskVal}
entries are passed to the psLib statistics function used to calculate
the ensemble statistics.  This function must be defined for the
following types: \code{psS8}, \code{psU16}, \code{psF32},
\code{psF64}.

\begin{prototype}
long psImageCountPixelMask (psImage *mask, psRegion region, psMaskValue value);
\end{prototype}
This function returns the number of pixels in the image region which
satisfy any of the mask bits.  An error (eg, invalid image, invalid
region) results in a return value of -1.

\begin{prototype}
psPolynomial2D *psImageFitPolynomial(psPolynomial2D *coeffs, const psImage *input);
\end{prototype}
Fit a 2-D Chebychev polynomial surface to an image.  The input
structure \code{coeffs} contains the desired order and terms of
interest.  This function must be defined for the
following types: \code{psS8}, \code{psU16}, \code{psF32}, \code{psF64}.

\begin{prototype}
psImage *psImageEvalPolynomial(psImage *input, const psPolynomial2D *coeffs);
\end{prototype}
Evaluate a 2-D polynomial surface for the image pixels.  Given the
input polynomial coefficients, set the image pixel values on the basis
of the polynomial function.  This function must be defined for the
following types: \code{psS8}, \code{psU16}, \code{psF32}, \code{psF64}.

\begin{prototype}
complex psImagePixelInterpolate(const psImage *input, float x, float y,
                              const psImage *mask, psMaskType maskVal,
                              complex unexposedValue, psImageInterpolateMode mode);
\end{prototype}
Perform interpolation of image pixel values to the given fractional
coordinate \code{x,y}.  The function returns the interpolated value of
the image at the given fractional pixel coordinates, based on the
specified interpolation \code{mode}.  The \code{mask} allows pixels to
be excluded if their corresponding mask pixel value matches the value
of \code{maskVal}.  This function will likely be implemented as a
macro for processing speed.  It may also be necessary to define a
setup macro which pre-calculates certain values which would be reused
in a loop.

\subsubsection{Image Pixel Manipulations}

\begin{prototype}
int psImageClip(psImage *input, double min, double vmin, double max, double vmax);
\end{prototype}
Clip image values outside of range to given values.  All pixels with
values \code{< min} are set to the value \code{vmin}. All pixels with
values \code{> max} are set to the value \code{vmax}. Returns the
number of clipped pixels.  This function must be defined for the
following types: \code{psU8}, \code{psU16}, \code{psS8}, \code{psS16},
\code{psF32}, \code{psF64}, \code{psC32}, \code{psC64}.  The arguments
(\code{min}, \code{max}, etc) must be cast to the appropriate types to
match the image data.  If the input parameters \code{vmin} or
\code{vmax} are out of bounds for the image pixel type, the function
must raise an error.  It is not an error for \code{min} or \code{max}
to be out of range.  In the case of complex numbers, the input
parameters \code{min} and \code{max} must be compared against the
absolute value of the pixel values.  The values to which complex image
data are clipped employ the provided value for the real component and
0.0 for the imaginary component.

\begin{prototype}
int psImageClipComplexRegion(psImage *input, complex min, complex vmin,
                             complex max, complex vmax);
\end{prototype}
Clip image values outside of range to given values.  All pixels with
values \code{< min} are set to the value \code{vmin}. All pixels with
values \code{> max} are set to the value \code{vmax}. Returns the
number of clipped pixels.  This function must be defined for the
following types: \code{psC32}, \code{psC64}.  The arguments
(\code{min}, \code{max}, etc) define a rectangular region in complex
space; data values inside this regions are unchanged while those
outside are set to either \code{vmax} (if either their real or
imaginary portions are greater than the corresponding values of
\code{max}) or \code{vmin} (in all other cases).  If the input
parameters \code{vmin} or \code{vmax} are out of bounds for the image
pixel type, the function must raise an error.  It is not an error for
\code{min} or \code{max} to be out of range.

\begin{prototype}
int psImageClipNaN(psImage *input, float value);
\end{prototype}
Clip \code{NaN} image pixels to given value.  Pixels with \code{NaN},
\code{+Inf} or \code{-Inf} values are set to the specified value.
Returns the number of clipped pixels.  This function must be defined
for the following types: \code{psF32}, \code{psF64}, \code{psC32},
\code{psC64}.  In the case of complex values, if either the real or
imaginary part have the value \code{NaN}, then that component will be
set to the specified value.

\begin{prototype}
int psImageOverlaySection(psImage *image, const psImage *overlay, 
                          int x0, int y0, const char *op);
\end{prototype}
Overlay subregion of image with another image.  Replace the pixels in
the \code{image} which correspond to the pixels in \code{overlay} with
values derived from the values in \code{image} and \code{overlay}
based on the given operator \code{op}.  Valid operators are \code{=}
(set image value to overlay value), \code{+} (add overlay value to
image value), \code{-} (subtract overlay from image), \code{*}
(multiply overlay times image), \code{/} (divide image by overlay).
This function must be defined for the following types: \code{psU8},
\code{psU16}, \code{psS8}, \code{psS16}, \code{psF32}, \code{psF64},
\code{psC32}, \code{psC64}.  The two input images must have the same
datatype and the output image must be constructed with that same
datatype.  The return value shall be the number of pixels overlaid.

\subsubsection{Mask operations}

\begin{prototype}
psImage *psImageGrowMask(psImage *out, const psImage *in, psU16 maskVal,
                         unsigned int growSize, psU16 growVal);
\end{prototype}

\code{psImageGrowMask} grows specified values on the input mask image,
\code{in}, returning the result.  If \code{out} is \code{NULL}, then a
new image of the same type and dimension as \code{in} shall be
allocated and returned; otherwise \code{out} shall be modified.  If
\code{out} is non-\code{NULL} and does not have the same size and type
as \code{in}, the function shall generate an error and return
\code{NULL}.  Pixels in the \code{in} image within \code{growSize}
pixels (either horizontal or vertical) of a pixel which matches the
\code{maskVal} shall have the corresponding pixel in the \code{out}
image set to the \code{growValue}.  The function must be defined for
the following types: \code{psU8}, \code{psU16}.  Note that
\code{maskVal} and \code{growVal} are not of type \code{psMaskType}, but
\code{psU16} so that they match the maximum required type for the
\code{in} image.

\begin{prototype}
psImage *psPixelsToMask(psImage *out, const psPixels *pixels, psRegion region, psMaskType maskVal);
psPixels *psPixelsFromMask(psPixels *out, const psImage *mask, psMaskType maskVal);
\end{prototype}

\code{psPixelsToMask} shall return an image of type \code{psMaskType} with
the \code{pixels} lying within the specified \code{region} set to the
\code{maskVal}.  The \code{out} image shall be modified if supplied,
or allocated and returned if \code{NULL}.  The size of the output
image shall be \code{region.x1 - region.x0} by \code{region.y1 -
region.y0}, with \code{out->x0 = region.x0} and \code{out->y0 =
region.y0}.  In the event that either of \code{pixels} or
\code{region} are \code{NULL}, the function shall generate an error
and return \code{NULL}.

\code{psMaskToPixels} shall return a \code{psPixels} containing the
coordinates in the \code{mask} that match the \code{maskVal}.  The
\code{out} pixel list shall be modified if supplied, or allocated and
returned if \code{NULL}.  In the event that \code{mask} is
\code{NULL}, the function shall generate an error and return
\code{NULL}.

\begin{prototype}
void psImageMaskRegion(psImage *image, psRegion region, const char *op, psMaskType maskValue);
void psImageKeepRegion(psImage *image, psRegion region, const char *op, psMaskType maskValue);
\end{prototype}

These two complementary functions set bit specified bits
(\code{maskValue}) in the mask \code{image} interior to or exterior to
the specified \code{region}.  The first function,
\code{psImageMaskRegion}, sets the bits inside the region, ignoring
pixels outside.  The second, \code{psImageKeepRegion}, sets the bits
outside the region, ignoring pixels inside.  The pixel values are set
by combining the existing pixel value and the given \code{maskValue}
with a logical operation.  The allowed operations are \code{=},
\code{AND}, \code{OR} and \code{XOR}.

\begin{prototype}
void psImageMaskCircle(psImage *image, double x, double y, double radius, const char *op,
                       psMaskType maskValue);
void psImageKeepCircle(psImage *image, double x, double y, double radius, const char *op,
                       psMaskType maskValue);
\end{prototype}

These two complementary functions set bit specified bits
(\code{maskValue}) in the mask \code{image} interior to or exterior to
the specified circle, defined by the center coordinates \code{x,y} and
a \code{radius}.  The first function, \code{psImageMaskCircle}, sets
the bits inside the circle, ignoring pixels outside.  The second,
\code{psImageKeepCircle}, sets the bits outside the circle, ignoring
pixels inside.  The pixel values are set by combining the existing
pixel value and the given \code{maskValue} with a logical operation.
The allowed operations are \code{=}, \code{AND}, \code{OR} and
\code{XOR}.

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

\subsection{Vector and Image Arithmetic}
\label{sec:arithmetic}

We will need to be able to perform various operations on vectors and
images, e.g.\ dividing one image by another, subtracting a vector from
an image, etc.  Both binary operations and unary operations are
required.  To avoid the burden of memorizing a plethora of APIs, we
specify two generic APIs for the binary and unary operations.

\begin{prototype}
psMathType *psBinaryOp(psPtr out, const psPtr in1, const char *op, const psPtr in2);
psMathType *psUnaryOp(psPtr out, const psPtr in, const char *op);
\end{prototype}
These functions determine the type of the operands on the basis of
their \code{psMathType} elements, which always are the first elements.
Note that these functions return a pointer to the appropriate type for
the operation.  Since the result may be cast to \code{psMathType}, the
resulting type may be determined by examining the return value.  It is
expected that the implementation of these functions will employ
pre-processor macros to perform the onerous task of creating the
loops.  An attempt to perform an arithmetic operation on an object of
dimension \code{PS_DIMEN_OTHER} should produce an error.  Operations
between data structures with different types (e.g., \code{psS32} and
\code{psF32}) are not allowed and must raise an error (it is the
responsibility of calling functions to perform type conversions).
Operations between data structures with incompatible sizes are not
allowed.  However, operations between data elements of different rank
(scalar, vector, image) are allowed, and defined below.  These
functions are valid for all data types \code{psU8}, \code{psU16},
\code{psS8}, \code{psS16}, \code{psF32}, \code{psF64}, \code{psC32},
\code{psC64}.

Binary operations between an image and a vector have a potential
ambiguity --- do the vector elements correspond to the rows or the
columns?  For this reason, we define two vector types: a ``vector''
(\code{PS_DIMEN_VECTOR}), and a ``transposed vector''
(\code{PS_DIMEN_TRANSV}).  We specify that a ``vector'', when involved
in binary operations on an image, acts on all rows of the image, while
a ``transposed vector'' in the same context acts on all columns.
Vectors, when created, will be created as ``vectors'', but may be
converted to ``transposed vectors'' by setting
\code{vector->type.dimen = PS_DIMEN_TRANSV}.

It is further desirable to allow scalar values to be used within these
functions.  These functions may take a pointer to type
\code{psScalar}, which is freed by the function. This allows one to
write the following lines to take the sine of the square of all pixels
in an image:
\begin{verbatim}
psImage A,B;

B = psBinaryOp (NULL, A, "^", psScalarAlloc(2)); 
psUnaryOp(B, B, "sin");
\end{verbatim}

Note that the \code{psUnaryOp} is performed on \code{B} in-place.

The list of required operators for \code{psBinaryOp} are:
\begin{itemize}
\item Addition ($+$)
\item Subtraction ($-$)
\item Multiplication ($*$)
\item Division ($/$)
\item Power (\^)
\item Minimum (min)
\item Maximum (max)
\item Logical or ($|$; integer type only)
\item Logical and (\&; integer type only)
\end{itemize}

The list of required operators for \code{psUnaryOp} are:
\begin{itemize}
\item Absolute value (abs)
\item Exponent (exp)
\item Natural Log (ln)
\item Power of 10 (ten)
\item Log (log)
\item Sine (sin and dsin)
\item Cosine (cos and dcos)
\item Tangent (tan and dtan)
\item Arcsine (asin and dasin)
\item Arccosine (acos and dacos)
\item Arctan (atan and datan)
\end{itemize}
The trigonometric operators prefixed with ``d'' refer to the standard
trigonometric operators acting on input that is in decimal degrees.

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


\subsection{Matrix operations and linear algebra}

In addition to the ability to perform image arithmetic
(\S\ref{sec:arithmetic}), we also require the ability to perform basic
linear algebra on matrices, in order to solve equations.  We use
\code{psImage} as a matrix, since it has the correct form.  We define
the following basic matrix operations:
\begin{itemize}
\item $LU$ Decompose a matrix, and solve for $x$ in $Ax = b$;
\item Invert a matrix;
\item Calculate a matrix determinant;
\item Perform matrix addition, subtraction and multiplication;
\item Transpose a matrix;
\item Get eigenvectors for a matrix; and
\item Convert a matrix to a vector.
\end{itemize}
The corresponding APIs follow.

In the event of an error, the matrix functions shall return \code{NULL}.

\begin{prototype}
psImage *psMatrixLUD(psImage *out, psVector **perm, const psImage *in);
psVector *psMatrixLUSolve(psVector *out, const psImage *LU, const psVector *RHS, const psVector *perm);
\end{prototype}
The above functions decompose a matrix, \code{in}, into its $LU$
representation (\code{psMatrixLUD}, which returns the decomposed
matrix), and uses the decomposed matrix to solve for $x$ in the
equation $Ax = b$.  In this case, the $LU$ decomposed matrix $A$ is
specified as \code{LU}, and $b$ is \code{RHS}; the solution vector for
$x$ is returned.  These functions are specified for data types
\code{psF32, psF64}.  The output and input vectors and images must all
have the same data type.

The GSL routines require the use of a permutation vector, \code{perm}.
This vector shall be created by \code{psMatrixLUD}, and the user need
only pass this to \code{psMatrixLUSolve} before destroying it in the
standard manner.  In order to avoid memory leaks, \code{perm} must be
\code{NULL} on calling \code{psMatrixLUD}.

\begin{prototype}
psImage *psMatrixInvert(psImage *out, const psImage *in, float *determinant);
\end{prototype}
\code{psMatrixInvert} returns the inverse of the specified matrix
(\code{in}), along with the \code{determinant}, if the \code{float}
pointer is non-NULL.  The input and output images must have the same
type.  This function is specified for data types \code{psF32, psF64}.

The matrix determinant may be calculated using
\code{psMatrixDeterminant}, which simply returns the determinant for
the specified matrix, \code{in}.  This function is specified for data
types \code{psF32, psF64}.
\begin{prototype}
float psMatrixDeterminant(const psImage *in);
\end{prototype}

Matrix multiplication is supported through \code{psMatrixMultiply}.  This function is
specified for data types \code{psF32, psF64}.
\begin{prototype}
psImage *psMatrixMultiply(psImage *out, const psImage *in1, const psImage *in2);
\end{prototype}

The transpose of an input matrix, \code{in}, is returned by
\code{psMatrixTranspose}, optionally into the provided matrix
\code{out}.  This function is specified for data types \code{psF32,psF64}.
\begin{prototype}
psImage *psMatrixTranspose(psImage *out, const psImage *in);
\end{prototype}

Eigenvectors of a matrix are calculated by
\code{psMatrixEigenvectors}.  This function is specified for data
types \code{psF32, psF64}.  Input and output data types should match.
\begin{prototype}
psImage *psMatrixEigenvectors(psImage *out, const psImage *in);
\end{prototype}
\tbd{Should this return an array of vectors?  Specified here as
currently implemented by MHPCC.}

Finally, we specify two functions to convert between matrices and
vectors.  This allows the use of vectors in matrix arithmetic, since a
vector can be converted to a matrix, utilized, and the result can be
converted back to a vector if required.  This function is specified
for data types \code{psF32, psF64}.  Input and output data types
should match.
\begin{prototype}
psVector *psMatrixToVector(psVector *outVector, const psImage *inImage);
psImage *psVectorToMatrix(psImage *outImage, const psVector *inVector);
\end{prototype}

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

\subsection{(Fast) Fourier Transforms}

We require the ability to calculate the (fast) Fourier transforms of
floating-point one-dimensional vectors and two-dimensional images.  We
expect that these will be implemented through wrapping an external
library.  We define the following APIs to support FFT operations on vectors:

\begin{prototype}
psVector *psVectorFFT(psVector *out, const psVector *in, psFFTFlags direction);
psVector *psVectorReal(psVector *out, const psVector *in);
psVector *psVectorImaginary(psVector *out, const psVector *in);
psVector *psVectorComplex(psVector *out, const psVector *real, const psVector *imag);
psVector *psVectorConjugate(psVector *out, const psVector *in);
psVector *psVectorPowerSpectrum(psVector *out, const psVector *in);
\end{prototype}

The forward and reverse FFT is calculated using \code{psVectorFFT},
which takes as input the vector of interest (\code{in}) and the
direction (\code{direction}), which is specified by an enumerated type
defined below.  The input vector may be of type \code{psF32} or
\code{psC32}, the result is always \code{psC32}.  If the input vector
is \code{psF32}, the direction must be forward.  Neither the forward
or inverse transforms must multiply by $1/N$ (or $1/N^{1/2}$), and so
it falls to the responsibility of the user to multiply a vector that
has been forward- and reverse-transformed by $1/N$.

Conversions between complex and real vectors requires the functions
\code{psVectorReal}, which returns the real part (\code{psF32}) of the
complex vector \code{in}, \code{psVectorImaginary}, which returns the
the magnitude of the imaginary part (\code{psF32}) of the complex
vector \code{in}, and \code{psVectorComplex}, which constructs a
complex vector (\code{psC32}) from the real and imaginary components
\code{real} and \code{imag}, both of type \code{psF32}.

We also specify the additional utility functions
\code{psVectorConjugate} and \code{psVectorPowerSpectrum} which
construct the complex conjugate (\code{psC32}) and the magnitude
(\code{psF32}) vectors of the input complex vector (\code{psC32}).

The direction of an FFT performed by \code{psVectorFFT} is specified
by an enumerated type, \code{psFFTFlags}:

\begin{datatype}
/** Specify direction of FFT, and if the result is real or not */
typedef enum {
    PS_FFT_FORWARD = 1,                 ///< psImageFFT/psVectorFFT should perform a forward FFT.
    PS_FFT_REVERSE = 2,                 ///< psImageFFT/psVectorFFT should perform a reverse FFT.
    PS_FFT_REAL_RESULT = 4              ///< psImageFFT/psVectorFFT should return a real image. This is valid
                                        ///< for only reverse FFT, i.e., the psImageFFT/psVectorFFT flag
                                        ///< parameter should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT.
} psFFTFlags;
\end{datatype}

The entry \code{PS_FFT_REAL_RESULT} means that the output of an
inverse FFT is to be purely real.  An example of its use is:
\begin{verbatim}
out = psImageFFT(in, PS_FFT_REVERSE | PS_FFT_REAL_RESULT);
\end{verbatim}

The output from a FFT must be of the same size as the input, and the
size of the power spectrum equal to $N/2 + 1$, where $N$ is the number
of elements in the input.  In the future, if this adversely affects
performance, this will be revised so that the redundant information
will be neglected.

In analogy with the vector FFT operations, we also define matching
image operations as follows:
\begin{prototype}
psImage *psImageFFT(psImage *out, const psImage *image, psFFTFlags direction);
psImage *psImageReal(psImage *out, const psImage *in);
psImage *psImageImaginary(psImage *out, const psImage *in);
psImage *psImageComplex(psImage *out, const psImage *real, const psImage *imag);
psImage *psImageConjugate(psImage *out, const psImage *in);
psImage *psImagePowerSpectrum(psImage *out, const psImage *in);
\end{prototype}

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

\subsection{Convolution}

Convolution will be an essential operation for the IPP.  For example,
if images on the sky are obtained in Orthogonal Transfer (OT) mode,
then the calibration frames used to correct them must be convolved by
a kernel derived from the list of OT shifts made during the exposure.
Also, convolution is also required for object detection and
PSF-matching.

\subsubsection{Basic Image Smoothing}

\begin{prototype}
void psImageSmooth (psImage *image, double sigma, double Nsigma);
\end{prototype}

The quickest way to smooth an image is to smooth by parts, using 1D
gaussian smoothing in X and Y independently.  \code{psImageSmooth}
applies a single parameter Gaussian (circularly symmetric) by
smoothing first in the X and then in the Y directions with just a
vector.  In general, for a Gaussian of dimension N, this process is 2N
faster than for the 2D convolution defined below.  The arguments to
this function include the image to be smoothed, the width of the
smoothing kernel in pixels (\code{sigma}), and the size of the
smoothing box in sigmas.  

\subsubsection{Kernel definition}

In order to perform a convolution, we need to define the convolution
kernel.  We need a more general object than a \code{psImage} so that
we can incorporate the offset from the $(0,0)$ pixel to the $(0,0)$
value of the kernel.  It might be convenient to allow both positive
and negative indices to convey the positive and negative shifts.  One
might consider setting the \code{x0} and \code{y0} members of a
\code{psImage} to the appropriate offsets, but this is not the purpose
of these members, and doing so may affect the behavior of other
\code{psImage} operations.  We define a \code{psKernel}:

\begin{datatype}
/** A convolution kernel */
typedef struct {
    psImage *image;                     ///< Kernel data, in the form of an image
    int xMin;                           ///< Most negative indices
    int yMin;                           ///< Most negative indices
    int xMax;                           ///< Most positive indices
    int yMax;                           ///< Most positive indices
    float **kernel;                     ///< Pointer to the kernel data
    float **p_kernelRows;               ///< Pointer to the rows of the kernel data
} psKernel;
\end{datatype}

The kernel data is carried primarily by the \code{data} member which
is a normal \code{psImage}.  In order to allow negative indices, we
add two additional members.  \code{kernelRows} is an array of pointers
to \code{float}; these pointers point into the \code{psImage} data,
offset by the desired column offset.  \code{kernel} point into the
\code{kernelRows}, offset by the desired row offset.  This
construction allows the \code{kernel} member to use negative indices,
while preserving the location of \code{psMemBlock}s relative to
allocated memory.

The maximum extent of the kernel shifts shall be defined by the
\code{xMin}, \code{xMax}, \code{yMin} and \code{yMax} members.  Note
that \code{xMin} and \code{yMin}, under normal circumstances, should
be negative numbers.  That is, \code{myKernel->kernel[-3][-2]} may be
defined if \code{yMin} and \code{xMin} are equal to or more negative
than -3 and -2, respectively.

Of course, we require the appropriate constructor:
\begin{prototype}
psKernel *psKernelAlloc(int xMin, int xMax, int yMin, int yMax);
\end{prototype}

\code{psKernelAlloc} shall allocate a kernel.  In the event that one
of the minimum values is greater than the corresponding maximum value,
the function shall generate a warning, and the offending values shall
be exchanged.

\subsubsection{Generation of a convolution kernel}

Given a list of values (e.g., shifts made in the course of OT
guiding), \code{psKernelGenerate} shall return the appropriate kernel.
The API shall be the following:
\begin{prototype}
psKernel *psKernelGenerate(const psVector *tShifts, const psVector *xShifts,
                           const psVector *yShifts, bool relative);
\end{prototype}

The vectors \code{xShifts} and \code{yShifts}, which are a list of
shifts made at the times \code{tShifts}, are used to construct the
appropriate kernel.  If \code{relative} is \code{true}, then each
shift is to be interpreted relative the shift made before; if
\code{relative} is \code{false}, then the shifts are to be interpreted
relative to some starting point.  The elements of the vectors should
be of an integer type; otherwise the values shall be truncated to
integers.  The output kernel shall be normalized such that the sum
over the kernel is unity.

If the vectors are not all of the same number of elements, then the
function shall generate an error and return \code{NULL}.

\subsubsection{Convolve an image with the kernel}

Given an input image and the convolution kernel,
\code{psImageConvolve} shall convolve the input image,
\code{in}, with the kernel, \code{kernel} and return the convolved
image, \code{out}.  The API shall be the following:
\begin{prototype}
psImage *psImageConvolve(psImage *out, const psImage *in, const psKernel *kernel, bool direct);
\end{prototype}

Two methods shall be available for the convolution: if \code{direct}
is \code{true}, then the convolution shall be performed in real space
(appropriate for small kernels); otherwise, the convolution shall be
performed using Fast Fourier Transforms (FFTs; appropriate for larger
kernels).  The latter option involves padding the input image, copying
the kernel into an image of the same size as the padded input image,
performing an FFT on each, multiplying the FFTs, and performing an
inverse FFT before trimming the image back to the original size.

In the event that \code{out} is \code{NULL}, a new \code{psImage}
shall be allocated and returned.

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

\subsection{Random Numbers}

Many applications involve random numbers, with the particular
distribution depending upon the application.  We will wrap the GSL
implementation.

We define a \code{psRandom} type, which will allow us to carry around
pre-computed random numbers, if required.  For the time being, we only
specify a single random number generator (RNG) in \code{psRandomType}
(that provided by \code{gsl_rng_taus}), but we leave this open to
expansion in the future, depending upon user requirements.

\begin{datatype}
typedef enum {
    PS_RANDOM_TAUS                      ///< A maximally equidistributed combined Tausworthe generator
} psRandomType;

typedef struct {
    psRandomType type;                  ///< The type of RNG
    gsl_rng *gsl;                       ///< The RNG itself
} psRandom;
\end{datatype}

We require the ability to seed the random number generator (RNG) with
a known value, so that bugs in modules which rely upon random numbers
may be reproduced.

\begin{prototype}
psRandom *psRandomAlloc(psRandomType type, unsigned long seed);
void psRandomReset(psRandom *rand, unsigned long seed);
\end{prototype}

\code{psRandomAlloc} shall construct a new instance of \code{psRandom}
of the given \code{type}, and seed it with the given \code{seed}.  The
\code{seed} is specified as an \code{unsigned long} so that the system
clock may be used to set it.  If the \code{seed} is zero, then the RNG
shall be seeded from \code{/dev/urandom} if it exists, or otherwise
from the system clock, and the particular seed shall be printed using
\code{psLogMsg} with a level of \code{PS_LOG_INFO}.

\begin{prototype}
double psRandomUniform(const psRandom *r);
double psRandomGaussian(const psRandom *r);
double psRandomPoisson(const psRandom *r, double mean);
\end{prototype}

\code{psRandomUniform} shall return random numbers uniformly
distributed on $[0,1)$, using \code{gsl_rng_uniform}.
\code{psRandomGaussian} shall return random numbers distributed on a
Gaussian deviate, $N(0,1)$, using \code{gsl_ran_gaussian}.
\code{psRandomPoisson} shall return random numbers distributed on a
Poisson distribution with the given \code{mean} using
\code{gsl_ran_poisson}.

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

\pagebreak 
\section{Astronomy-Related Functions}

The previous sections of this document defined basic functionality
needed by a wide range of scientific programming.  The rest of this
document concentrates on aspects of the library which are specific to
astronomy.  Some basic, relatively simple astronomy-specific functions
are required which will serve as the foundation for building the
higher level modules.  These functions are not expected to cover every
foreseeable function, but will serve as the building blocks of more
complicated processing functions.  

We require functions covering each of the following areas:
\begin{itemize}
\item Dates and times
\item Detector and sky positions
\item Astrometry
\item Photometry
\item Astronomical objects
\end{itemize}

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

\subsection{Dates and times}
\label{sec:timespec}

\subsubsection{Overview}

We require a collection of functions to manipulate time data.  These
operations primarily consist of conversions between specific time
formats.  Internally, PSLib handles times as a structure similar to
the POSIX \code{struct timeval} struct which has been extended to
track the time system being represented.

\subsubsection{Initialization and Finalization}

The conversions between various time standards, and the polar motion
requires importing external data (``time tables''), the locations of
which are specified in a configuration file.

\begin{prototype}
void psTimeInitialize(const char *timeConfig);
\end{prototype}

\code{psTimeInitialize} shall read the \code{psTime} configuration
file (\code{timeConfig}) and set up the appropriate
\code{psTimeTable}s and predictions.

{\bf Calls to \code{psTime} functions that require the time tables
before calling \code{psTimeInitialize} shall result in an error.}

\begin{prototype}
void psTimeFinalize(void);
\end{prototype}

\code{psTimeFinalize} shall free memory that was allocated by
\code{psTime} functions, allowing a subsequent search for leaked
memory (even memory marked as \code{persistent}).

\subsubsection{Data Types}

\begin{datatype}
typedef enum {
    PS_TIME_TAI,      ///< seconds since 1970-01-01T00:00:00Z (Gregorian), SI seconds
    PS_TIME_UTC,      ///< seconds since 1970-01-01T00:00:00Z (Gregorian), SI seconds + leapseconds
    PS_TIME_UT1,      ///< seconds since 1970-01-01T00:00:00Z (Gregorian), variable length seconds
    PS_TIME_TT,       ///< seconds since 1970-01-01T00:00:00Z (Gregorian), SI seconds
} psTimeType;

typedef enum {
    PS_IERS_A,                      ///< IERS Bulletin A
    PS_IERS_B,                      ///< IERS Bulletin B
} psTimeBulletin;

typedef struct {
    psS64            sec;           ///< seconds, negative values represent dates before 1970
    psU32            nsec;          ///< nanoseconds
    bool             leapsecond;    ///< if time falls on a UTC leapsecond
    psTimeType       type;          ///< type of time
} psTime;
\end{datatype}

\subsubsection{Constructors}

To allocate a new, initialized to zero, \code{psTime}:

\begin{prototype}
psTime *psTimeAlloc(psTimeType type);
\end{prototype}

To allocate a new \code{psTime} set to the current time (in the given system):

\begin{prototype}
psTime *psTimeGetNow(psTimeType type);
\end{prototype}

\subsubsection{Time Conversion}

Converting between the \code{psTime} time systems is done with:

\begin{prototype}
psTime *psTimeConvert(psTime *time, psTimeType type);
\end{prototype}

This function may be used to convert between the various \code{psTimeType} time
representations.  The \code{time} is modified and returned.  Conversion between
all of the \emph{SI} length second systems should be implemented by first
converting to TAI and then to the destination system.  UT1 is a special case
for conversion as it uses variable length seconds.  Conversation to UT1, via
TAI, is allowed but conversion \emph{from} UT1 is currently forbidden.

The following function converts to Local Mean Sidereal Time.  It is
necessary to provide the local longitude (specified in radians,
positive East of Greenwich) as well:

\begin{prototype}
double psTimeToLMST(psTime *time, double longitude);
\end{prototype}

The function may accept any of the \code{psTimeType} types with
\emph{SI} length seconds.  The \code{time} is modified (rationlized to
TAI units) and returned.  Note that this function requires the value
$UT1-UTC$ (see \code{psTimeGetUT1Delta()}) and should use the UT1
values interpolated from IERS bulletin B.

The following utility function encapsulates the PSLib mechanism to
extract the value of $UT1-UTC$ from the IERS Time Tables:

\begin{prototype}
double psTimeGetUT1Delta(const psTime *time, psTimeBulletin bulletin);
\end{prototype}

The following function provides tidal corrections to UT1-UTC, using
the Ray model of Simon et al (REF).
\begin{prototype}
psTime *psTime_TideUT1Corr(const psTime *time);
\end{prototype}

Leap seconds are added to UTC in order to keep it within $0.9s$ of UT1
(which is defined relative to the Earth's rotation, and hence is
useful for astronomical purposes).  The following function calculates
the absolute number of leap seconds different between two times.

\begin{prototype}
long psTimeLeapSecondDelta(const psTime *time1, const psTime *time2);
\end{prototype}

The following function accepts \code{PS_TIME_UTC} objects and
determines if the time is potentially a leapsecond, returning
\code{TRUE} if so.

\begin{prototype}
bool psTimeIsLeapSecond(const psTime *utc);
\end{prototype}

\subsubsection{External Date and Time Formats}

A collection of functions convert from the \code{psTime} types to various
external formats.  The ISO8601 format\footnote{ISO8601:2000(E) -
http://www.pvv.ntnu.no/~nsaa/8601v2000.pdf} we will be using is
"YYYY-MM-DDThh:mm:ss,sZ".  \emph{Note that the ISO8601 format is only valid for
the years $0000 .. 9999$.  Values out side that range should cause a
\code{NULL} pointer to be returned.}

\begin{prototype}
double psTimeToJD(const psTime *time);
double psTimeToMJD(const psTime *time);
psString psTimeToISO(const psTime *time);
struct timeval *psTimeToTimeval(const psTime *time);
\end{prototype}

The \code{psTimeToISO()} function converts \code{PS_TIME_UTC} objects
with the \code{leapsecond} flag set to represent the number of seconds
as ``:60''.

\subsubsection{Date and Time Parsing}

A related collection of functions convert from external formats to a
\code{psTime} type.

\begin{prototype}
psTime *psTimeFromJD(double jd);
psTime *psTimeFromMJD(double mjd);
psTime *psTimeFromISO(const char *input, psTimeType type);
psTime *psTimeFromTimeval(const struct timeval *input);
psTime *psTimeFromUTC(psS64 sec, psU32 nsec, bool leapsecond);
psTime *psTimeFromTT(psS64 sec, psU32 nsec);
\end{prototype}

Where \code{psTimeFromUTC()} will validate that it is possible for the input
time to land on a UTC leapsecond (using \code{psTimeIsLeapSecond()}) if
\code{leapsecond} is set to \code{true}.

\subsubsection{Date and Time Math}

\begin{prototype}
psTime *psTimeMath(const psTime *time, double delta);
double psTimeDelta(const psTime *time1, const psTime *time2);
\end{prototype}

\code{psTimeMath} adds the \code{delta} (in seconds; may be negative) to a
\code{psTime}.  \code{PS_TIME_UTC} objects will be converted to
\code{PS_TIME_TAI}, before applying the delta, and then back to
\code{PS_TIME_UTC}.  \code{psTimeDelta} gives the difference between times
\code{time1} and \code{time2} (which may be negative).

The API documentation should note that when handling UT1, date math is allowed
on the UT1 system but that the seconds are not of \emph{SI} length (this is why
conversion from UT1 is correctly forbidden).  It should also be noted that with
\code{psTimeDelta()} it is possible to overflow the dynamic range of a
\code{psS64} (the type of \code{psTime.sec}).

Note that in both these functions the difference between two times is
the elapsed number of seconds, including leap seconds.  For example,
if we add 30 seconds to 1998-12-31T23:59:45Z, the result is
1999-01-01T00:00:14Z, since a leap second was introduced at
1999-01-01T00:00:00Z.

Time math may only be done on \code{psTime} objects of the same type,
and except in the case of UT1, the functions shall internally convert
the \code{psTime} inputs to TAI before performing the math; this
ensures that leap seconds are correctly counted.

The type of the time returned by \code{psTimeMath} shall be the same
as that of the input \code{time}.

\subsubsection{Time Tables}

The offset of UTC from UT1, $\Delta UT = UT1-UTC$, as well as the polar motion,
$x_p$ and $y_p$, may be determined from table lookups.  Tables are available
covering different time periods and with different time resolution, and so it
is important to be able to utilize multiple tables.  Some tables may be found
at:

\begin{itemize}
\item IERS Bulletin A \& B (1 year ago $\rightarrow$ now + $\sim$3 months)
\begin{itemize}
\item \code{ftp://maia.usno.navy.mil/ser7/finals.daily}
\end{itemize}

\item IERS Bulletin A \& B (1973-1-2 $\rightarrow$ now + $\sim$1 year)
\begin{itemize}
\item \code{ftp://maia.usno.navy.mil/ser7/finals.all}
\end{itemize}

\item $TAI - UTC$
\begin{itemize}
\item \code{ftp://maia.usno.navy.mil/ser7/tai-utc.dat}
\end{itemize}

\item $TAI -UTC$ (contains estimates prior to 1972)
\begin{itemize}
\item \code{http://hpiers.obspm.fr/eoppc/eop/eopc01/eopc01.1900-2004}
\end{itemize}
\end{itemize}

The format of \code{finals.daily} and \code{finals.all} is documented at
\code{ftp://maia.usno.navy.mil/ser7/readme.finals}.

The tables shall reside on local disk in known locations (i.e., there is no
need that they are downloaded from the internet and parsed by PSLib).  The
format of these files shall be simple, for speed in reading.  Each line shall
contain the date in MJD and the following values from both the A \& B IERS
bulletins: $x_p$ (in arcseconds), $y_p$ (in arcseconds) and $\Delta$UT (in
seconds).  This format must be readable by \code{psLookupTableRead}.  For
example:

\begin{verbatim}
# finals.daily, 2005-03-17
# ftp://maia.usno.navy.mil/ser7/finals.daily
# MJD    PM-IP  PM-X"   PM-Y"   UT1-YP  UT1-UTCs    PM-X"   PM-Y"   UT1-UTCs
#        A      A       A       A       A           B       B       B
53403.00 I      .089198 .207785 I       -.5218847   .089220 .207360 -.5219040
53404.00 I      .086549 .206873 I       -.5224164   .086670 .206850 -.5224260
53405.00 I      .083589 .206143 I       -.5227681
53406.00 I      .081541 .205865 I       -.5229582
\end{verbatim}

The location of these files, their priority order, and the ``from''
and ``to'' dates of applicability will be specified through metadata
(\S\ref{sec:timeMetadata}).

{\bf Calls to \code{psTime} functions that require the time tables
before calling \code{psTimeInitialize} shall result in an error.}

When a value is required, the tables shall shall be checked in
priority order to see if the date is within the range of applicability
for the table.  If a table is found that is applicable, then the
appropriate value shall be derived from linear interpolation between
the nearest entries in the table.  If no table is found that is
applicable, and the required date is later than those covered in the
tables, a warning shall be generated, and a pre-determined formula
shall be applied (\S\ref{sec:timeMetadata}).  For dates prior to those
covered in the tables, the function shall generate a warning and use
pre-determined values (\S\ref{sec:timeMetadata}).

\paragraph{Time Metadata}
\label{sec:timeMetadata}

The following metadata keys will be used in time calculations:

\begin{tabular}{l|l} \hline
Metadata key & Purpose \\ \hline
psLib.time.tables.dir   & Time table directory                                             \\
psLib.time.tables.files & Time table file names (space-delimited)                          \\
psLib.time.tables.from  & Time tables are valid from these MJDs (vector)                   \\
psLib.time.tables.to    & Time tables are valid to these MJDs (vector)                     \\
psLib.time.before.xp    & Value of XP for before the earliest MJD                          \\
psLib.time.before.yp    & Value of YP for before the earliest MJD                          \\
psLib.time.before.dut   & Value of UT1-UTC for before the earliest MJD                     \\
pslib.time.predict.xp   & A vector containing the $x_p$ prediction formula coefficients    \\
pslib.time.predict.yp   & A vector containing the $y_p$ prediction formula coefficients    \\
pslib.time.predict.mjd  & A value containing the MJD offset for temporary variables        \\
pslib.time.predict.dut  & A vector containing the UT1-UTC prediction formula coefficients  \\
\hline
\end{tabular}

These metadata keys shall reside in a configuration file
(\S\ref{sec:configspec}), \code{psTime.config}, which shall be loaded
into a \code{psMetadata} structure private to \code{psTime}, as part
of \code{psLibInit}.  An example of \code{psTime.config} follows:

\begin{verbatim}
# This configuration file specifies values required for time calculations by psLib.
psLib.time.tables.dir   STR     /home/panstarrs/psLib/config/                   # Directory for time tables
psLib.time.tables.n     U8      2                                               # Number of time tables
psLib.time.tables.files STR     bulletinA_09Sep2004.dat eopc01_1900_2004.40.dat # These are the file names
@psLib.time.tables.from F64     53258.0, 15020.0                                # Valid from these MJDs
@psLib.time.tables.to   F64     53622.0, 53258.0                                # Valid to these MJDs
psLib.time.before.xp    F64     0.0                     # Value of XP for before the earliest MJD
psLib.time.before.yp    F64     0.0                     # Value of YP for before the earliest MJD
psLib.time.before.dut   F64     0.0                     # Value of UT1-UTC for before the earliest MJD
                                                        
# Now follows formulae for predicting ahead of the most recent available table entry.
# xp = [0] + [1]*cos A + [2]*sin A + [3]*cos C + [4]*sin C
# yp = [0] + [1]*cos A + [2]*sin A + [3]*cos C + [4]*sin C
# A = 2*pi*(MJD - pslib.time.predict.mjd)/365.25
# C = 2*pi*(MJD - pslib.time.predict.mjd)/435.0
# dut = ut1-utc = [0] + [1]*(MJD - [2]) - (UT2-UT1)
# ut2-ut1 = 0.022 sin(2*pi*T) - 0.012 cos(2*pi*T) - 0.006 sin(4*pi*T) + 0.007 cos(4*pi*T)
# T = 2000.0 + (MJD - 51544.03)/365.2422
@psLib.time.predict.xp  F64     0.0569, 0.0555, -0.0200, 0.0513, 0.1483
@psLib.time.predict.yp  F64     0.3498, -0.0176, -0.0498, 0.1483, -0.0513
psLib.time.predict.mjd  F64     53257.0
@psLib.time.predict.dut F64     -0.4944, -0.00023, 53262.0
\end{verbatim}

A series of test inputs is contained in \S\ref{sec:timetest}.

\subsection{Timers}

It is useful to be able to time an operation.  For this purpose, we specify
the following:

\begin{prototype}
bool psTimerStart(char *name);
psF64 psTimerMark(char *name);
psF64 psTimerClear(char *name);
psF64 psTimerStop(void);
\end{prototype}

\code{psTimerStart} shall store the current time in a \code{psHash} of
timers, under the supplied \code{name}.  \code{psTimerMark} shall
return the elapsed time, in seconds, for the timer specified by
\code{name}.  \code{psTimerClear} resets the named timer.
\code{psTimerStop} shall free all memory associated with all timers,
so that no memory is leaked, and return the maximum time expended.

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

\subsection{Linear and Spherical Coordinates}

Both detector and sky positions will be used extensively in the IPP.
The first are linear coordinates which conform to Euclidean geometry
while the second are angular coordinates which define a position on
the sphere of the sky.  We put these into two structures,
\code{psPlane} and \code{psSphere}, respectively.  Partitioning these
two will enable error-checking.  An alternative representation for
angular positions is the 3-D unit vector.  These are used in
particular as part of spherical rotation calculations.  We define
\code{psCube} to represent such an element.
%
\begin{datatype}
typedef struct {
    double x;                           ///< x position
    double y;                           ///< y position
    double xErr;                        ///< Error in x position
    double yErr;                        ///< Error in y position
} psPlane;

typedef struct {
    double r;                           ///< RA
    double d;                           ///< Dec
    double rErr;                        ///< Error in RA
    double dErr;                        ///< Error in Dec
} psSphere;

typedef struct {
    double x;                           ///< cos (DEC) cos (RA) 
    double y;                           ///< cos (DEC) sin (RA) 
    double z;                           ///< sin (DEC)
    double xErr;                        ///< Error in x
    double yErr;                        ///< Error in y
    double zErr;                        ///< Error in z
} psCube;
\end{datatype}

Three major classes of coordinate transformations are necessary.
First, linear coordinates from one frame must be converted to linear
coordinates in a different frame of references.  Simple
transformations of this type are independent of other quantities of
the positions -- they are simply mapping between two linear spaces.
In practice, these transformations may often be a function of the
color or even magnitude of the imaged object.  The second type of
conversion is the transformation of linear coordinates to angular
coordinates and vice-versa.  This conversion depends on the desired
projection, and may represent the real mapping performed by the
telescope or may simply represent a convenient mechanism to display 3D
coordinates in useful forms.  The third conversion of interest is the
transformation of one set of spherical coordinates to another set.
Frequently in astronomy, these conversions consist only of rotations
between the two spherical coordinates systems, where the coordinates
of the pole and equatorial rotation between the two systems define the
transformation.  Conversions between standard coordinate systems such
as Galactic, Ecliptic, and various epochs of the Celestial coordinates
are represented by these spherical transformations.

Constructors for these are straight-forward:
\begin{prototype}
psPlane  *psPlaneAlloc(void);
psSphere *psSphereAlloc(void);
psCube   *psCubeAlloc(void);
\end{prototype}
Initialization of the structures is not necessary.

\subsubsection{Linear Coordinate Transformations}

We specify two types of transforms between coordinate systems.  The
first consists simply of two 2D polynomials to transform both
components -- the output coordinates depend only on the input
coordinates and no other quantities of objects at those coordinates.
The second consists of two 4D polynomials in which the output
coordinates are also specified to be a function of the color and
magnitude of the object with the given coordinates.  This type of
coordinate transformation is necessary to represent the
(color-dependent) optical distortions caused by the atmosphere and
camera optics, and the possibly effects of charge transfer
inefficiency.  We specify two structures to represent the coefficients
of these transformations:

\begin{datatype}
typedef struct {
    psPolynomial2D *x;
    psPolynomial2D *y;
} psPlaneTransform;
\end{datatype}

The \code{psPolynomial2D} structures represent polynomials of
arbitrary order as a function of two dimensions.  There is one of
these structures for each of the two output dimensions.  As an
example, consider the simple transformation from one linear coordinate
frame ($x,y$), e.g., on a CCD, to a second frame ($p,q$), e.g., on a
chip. If we have only first order terms in the transformation
\code{psPlaneTransform T}, the new coordinates would be represented by
the terms:
%
\begin{verbatim}
p = T.x->coeff[0][0] + x*T.x->coeff[1][0] + y*T.x->coeff[0][1];
q = T.y->coeff[0][0] + x*T.y->coeff[1][0] + y*T.y->coeff[0][1];
\end{verbatim}
%
where we have excluded the basic cross-term ($x \times y$) by using
the mask: \code{T.x->mask[1][1] = 1; T.y->mask[1][1] = 1;}

The \code{psPlaneDistort} represents an optical distortion.  The
lowest two terms are the $x$ and $y$ axis of the target system.  The
higher two terms may represent color and magnitude terms.
\begin{datatype}
typedef struct {
    psPolynomial4D *x;
    psPolynomial4D *y;
} psPlaneDistort;
\end{datatype}

Like \code{psPlaneTransform}, \code{psPlaneDistort} contains two
\code{psPolynomial4D} structures representing polynomials of
arbitrary order as a function of four, rather than two dimensions.
There is one of these structures for each of the two output
dimensions.  In this structure, the highest two dimensions could
represent a color and magnitude.  As an example, consider the simple
transformation from one linear coordinate frame ($x,y$), e.g., on a
CCD, of an object with color and magnitude ($c,m$) to a second frame
($p,q$), e.g., the focal plane. If we have only first order terms in
the transformation \code{psPlaneDistort T}, the new coordinates
would be represented by the terms:
%
\begin{verbatim}
p = T.x->coeff[0][0][0][0] + x*T.x->coeff[1][0][0][0] + y*T.x->coeff[0][1][0][0] 
  + c*T.x->coeff[0][0][1][0] + m*T.x->coeff[0][0][0][1]
q = T.y->coeff[0][0][0][0] + x*T.y->coeff[1][0][0][0] + y*T.y->coeff[0][1][0][0] 
  + c*T.y->coeff[0][0][1][0] + m*T.y->coeff[0][0][0][1]
\end{verbatim}
%
where we have again excluded the cross-term ($x \times y$) by using the
mask.

Each of \code{psPlaneTransform} and \code{psPlaneDistort} has an
appropriate allocator that takes the polynomial order in each
dimension.  Both the \code{x} and \code{y} polynomials shall be have
the same dimensions.
\begin{prototype}
psPlaneTransform *psPlaneTransformAlloc(int n1, int n2);
psPlaneDistort *psPlaneDistortAlloc(int n1, int n2, int n3, int n4);
\end{prototype}

We require corresponding functions to apply the transformations to a
specified coordinate \code{coords}:
%
\begin{prototype}
psPlane *psPlaneTransformApply(psPlane *out, 
                               const psPlaneTransform *transform, 
                               const psPlane *coords);
psPlane *psPlaneDistortApply(psPlane *out, 
                             const psPlaneDistort *distort, 
                             const psPlane *coords, 
                             float mag, float color);
\end{prototype}


The following functions perform operations on transformations:

\begin{prototype}
psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out, const psPlaneTransform *in,
                                         psRegion region, int nSamples);
psPlaneTransform *psPlaneTransformCombine(psPlaneTransform *out, const psPlaneTransform *trans1,
                                          const psPlaneTransform *trans2, psRegion region,
                                          int nSamples);
bool psPlaneTransformFit(psPlaneTransform *trans, const psArray *source, const psArray *dest,
                         int nRejIter, float sigmaClip);
\end{prototype}

\code{psPlaneTransformInvert} shall return the transformation that
inverts the given transformation.  It may assume that the input
transformation is one-to-one, and that the inverse transformation may
be specified through using polynomials of the same type and order as
the forward transformation.  In the event that the input
transformation is linear, an exact solution may be calculated;
otherwise \code{nSamples} samples in each axis, covering the region
specified by \code{region} shall be used as a grid to fit the best
inverse transformation.  The function shall return \code{NULL} if it
was unable to generate the inverse transformation; otherwise it shall
return the inverse transformation.  In the event that \code{out} is
\code{NULL}, a new \code{psPlaneTransform} shall be allocated and
returned.

\code{psPlaneTransformCombine} takes two transformations
(\code{trans1} and \code{trans2}) and returns a single transformation
that has the effect of performing \code{trans1} followed by
\code{trans2}.  In the event that the input transformation is linear,
an exact solution may be calculated; otherwise \code{nSamples} samples
in each axis, covering the region specified by \code{region} shall be
used as a grid to fit the best inverse transformation.  The function
shall return \code{NULL} if it was unable to generate the
transformation; otherwise it shall return the transformation.

\code{psPlaneTransformFit} takes two arrays containing matched
coordinates (i.e., coordinates in the \code{source} array correspond
to the coordinates in the \code{dest} array) and returns the
best-fitting transformation.  The \code{source} and \code{dest} will
contain \code{psCoord}s.  In the event that the number of coordinates
in each is not identical, the function shall generate a warning, and
extra coordinates in the longer of the two shall be ignored.  The
\code{trans} transform may not be \code{NULL}, since it specifies the
desired order, polynomial type and any polynomial terms to mask.
\code{nRejIter} rejection iterations shall be performed, wherein
coordinates lying more than \code{sigmaClip} standard deviations from
the fit shall be rejected.

\paragraph{Derivatives}

In order to simply calculate which pixels in one coordinate frame
overlap those in another coordinate frame tied by a
\code{psPlaneTransform}, we need to calculate the derivative of a
transformation with respect to each coordinate.

\begin{prototype}
psPlane *psPlaneTransformDeriv(psPlane *out, const psPlaneTransform *transformation, const psPlane *coord);
\end{prototype}

\code{psPlaneTransformDeriv} shall return the derivatives of the
\code{transformation} at the specified coordinates, \code{coord}.
Both the derivatives with respect to $x$ and $y$ shall be returned.
If \code{out} is non-\code{NULL}, it shall be modified and returned;
otherwise a new \code{psPlane} shall be allocated and returned.  In
the event that either \code{transformation} or \code{coord} are
\code{NULL}, the function shall generate an error and return
\code{NULL}.

\begin{prototype}
psPixels *psPixelsTransform(psPixels *out, const psPixels *input, const psPlaneTransform *inToOut);
\end{prototype}

\code{psPixelsTransform} shall generate a list of pixels in the output
coordinate frame that overlap the \code{input} pixels in the input
coordinate frame through the specified transformation, \code{inToOut}.
Note that this is more complicated than simply transforming the
\code{input} pixels, but requires the evaluation of the derivatives of
the transformation in order to obtain the list of all pixels in the
output coordinate frame that possibly overlap the pixel in the input
coordinate frame (assuming that $x' + \Delta x' = f(x) + \Delta x
\times \partial f(x) / \partial x$, where $x'$ is the coordinate in
the output coordinate frame, $\Delta x'$ is the size of the region in
the output coordinate frame in the positive direction that overlaps
the input pixel, $x$ is the coordinate in the input coordinate frame,
$f(x)$ is the transformation from the input to the output coordinate
frames, and $\Delta x$ is the size of a pixel; care should be taken
with half-pixel problems).  In the event that \code{input} or
\code{inToOut} are \code{NULL}, the function shall generate an error
and return \code{NULL}.

\subsubsection{Spherical Rotations}

Spherical rotations represent coordinate transformation in 3-D, as
well as the effects of precession and nutation.  We need spherical
rotations to convert between ICRS, Galactic and Ecliptic
coordinates, and to determine Alt-Az coordinates for sources.  All of
these basic spherical transformations represent rotations of the
spherical coordinate reference.  We specify a general transformation
function which takes a structure, \code{psSphereRot}, defining the
transformation between two spherical coordinate systems.  The
structure contains the elements of a quaternion to represent the
spherical rotational.  We define two allocators for
\code{psSphereRot}, one which defines the rotation in terms of the
coordinate of the pole and the rotation about that pole.  The other
defines the rotation from the elements of the quaternion.  We also
specify functions to manipulate \code{psSphereRot} in several useful
way.

\begin{datatype}
typedef struct {
    double q0;
    double q1;
    double q2;
    double q3;
} psSphereRot;
\end{datatype}

The constructor is defined as follows:
\begin{prototype}
psSphereRot *psSphereRotAlloc(double alphaP, double deltaP, double phiP);
\end{prototype}
where \code{alphaP} and \code{deltaP} define the coordinates in the
input system of the axis of rotation (the north pole of the output
system), while \code{phiP} defines the rotation about that pole.  This
last angle is also equal to 270\degree - $\phi_a$, where $\phi_a$ is
the longitude in the output system of the ascending node (equatorial
intersection between the two systems, e.g, the first point of Ares).

The \code{psSphereRot} may also be constructed by supplying the
elements of the quaternion to the following function:
\begin{prototype}
psSphereRot *psSphereRotQuat(double q0, double q1, double q2, double q3);
\end{prototype}
This function normalizes the quaternion, so the input elements need
not be normalized.

Spherical coordinates may be transformed by providing the
transformation and the coordinate in the input system to
\code{psSphereRot}.  The output pointer may be optionally supplied, or
if \code{NULL}, is allocated by the function.

\begin{prototype}
psSphere *psSphereRotApply(psSphere *out, const psSphereRot *transform, const psSphere *coord);
\end{prototype}

The following function combines two rotations, to produce a single
rotation which is the equivalent of applying the first rotation and
then the second.  The output rotation may be supplied, or will be
allocated if \code{NULL}.

\begin{prototype}
psSphereRot *psSphereRotCombine(psSphereRot *out, const psSphereRot *rot1, const psSphereRot *rot2)
\end{prototype}

The following function changes the given rotation to its inverse:

\begin{prototype}
psSphereRot *psSphereRotInvert(psSphereRot *rot)
\end{prototype}

The 3-vector representation of the angles (\code{psCube}) is needed to
implement these functions, and is useful in other circumstances as
well.  Two utility functions are provided to convert between the
angular and 3-vector representations:
\begin{prototype}
psSphere *psCubeToSphere(const psCube *cube);
psCube *psSphereToCube(const psSphere *sphere);
\end{prototype}

\subsubsection{Offsets}
We require a function to calculate the offset between two positions on
the sky, as well as a function to apply an offset to a position.  The
first determines the offset (RA,Dec) on the sky between two positions.
The second applies the given offset to the coordinate.  Both an offset
mode and an offset unit may be defined.  The mode may be either
\code{PS_SPHERICAL}, in which case the specified offset corresponds to
an offset in angles, or it may be \code{PS_LINEAR}, in which case the
offset corresponds to a linear offset in a local projection.  The
offset unit may be in one of \code{PS_ARCSEC}, \code{PS_ARCMIN},
\code{PS_DEGREE}, and \code{PS_RADIAN}, which specifies the units of
the offset only.

\begin{prototype}
psSphere *psSphereGetOffset(const psSphere *position1, 
                            const psSphere *position2, 
                            psSphereOffsetMode mode,
                            psSphereOffsetUnit unit);

psSphere *psSphereSetOffset(const psSphere *position,
                            const psSphere *offset,
                            psSphereOffsetMode mode,
                            psSphereOffsetUnit unit);
\end{prototype}

\begin{datatype}
typedef enum {
    PS_SPHERICAL,                       ///< Offset on a sphere
    PS_LINEAR                           ///< Linear offset
} psSphereOffsetMode;

typedef enum {
    PS_ARCSEC,                          ///< Arcseconds
    PS_ARCMIN,                          ///< Arcminutes
    PS_DEGREE,                          ///< Degrees
    PS_RADIAN                           ///< Radians
} psSphereOffsetUnit;
\end{datatype}
Note that these should propagate the errors appropriately.

\subsection{Celestial Coordinate Systems}

The following functions simply return the appropriate
\code{psSphereRot} to convert between predefined spherical
coordinate systems (i.e., ICRS, Ecliptic and Galactic).  These are
constructors as well as the above \code{psSphereRotAlloc}.
%
\begin{prototype}
psSphereRot *psSphereRotICRSToEcliptic(const psTime *time);
psSphereRot *psSphereRotEclipticToICRS(const psTime *time);
psSphereRot *psSphereRotICRSToGalactic(void);
psSphereRot *psSphereRotGalacticToICRS(void);
\end{prototype}

\subsection{Earth Orientation Calculations}

One of the critical sets of calculations in astronomy is the sequence
of steps needed to convert between the celestial coordinates of an
object and the observed coordinates of the object.  This problem is
best divided into two major components: transformation between the
celestial sphere and coordinates relative to the surface of the solid
earth, excluding the effects of the atmosphere, and compensation for
the effects of the atmosphere.  In this section, we address the first
of these two transformations: the Earth Orientation Calculations.

The Earth Orientation Calculations are further subdivided into several
steps, illustrated in Figure~\ref{CoordinateSystems} .  Celestial
coordinates are defined in the International Celestial Reference
System (ICRS), which has the solar barycenter as its reference
position and velocity.  The next coordinate system is the Geocentric
Celestial Reference System (GCRS), which uses the earth barycenter as
a reference.  The transformation between these two includes the
aberration due to the Earth's velocity, the parallax of the object,
which depends on both the Earth's position and the distance to the
object of interest, and the general relativistic correction for the
bending of light as it approaches the Earth.

The next set of transformations compensate for the 3-D rotation of
the Earth on various timescales, including the effects of precession,
nutation, and simple solid-body rotation.  These calculations can be
performed using different amounts of information for higher levels of
precision.  Since the Earth's rotation is constantly affected by
stochastic processes (weather, earthquakes, etc), these conversions
are constantly modified by observations reported by authoritative
sources such as the US Naval Observatory.  The target of this
transformation is the International Terrestrial Reference System
(ITRS), which is fixed with respect to the Earth's crust.  This
transformation is subdivided into slow precession and nutation
(yielding the coordinate system CIP/CEO), followed by the Earth's
rotation (yielding the coordinate system CIP/TEO), and finally
corrections for the short-period motion of the Earth's pole.  

\begin{figure}
\begin{center}
\psfig{file=earthrot.ps,height=4in}
\caption{Earth Orientation Coordinate Frames\label{CoordinateSystems}}
\end{center}
\end{figure}

\subsubsection{Transformation from ICRS to GCRS}

\tbd{we need a function to construct the direction and speed elements
  given the time}.

\tbd{supply the velocity as an un-normalized 3 vector?}

\tbd{MHPCC: please code this section as currently specified.  We will
  define a function, and algorithm, to return the current velocity
  vector given a time and position, which can be fed to this
  function}.

\paragraph{Aberration}
The following function calculates the \code{apparent} position of a
star, given its \code{actual} position and the velocity vector of the
observer, represented as a speed and a direction:
\begin{prototype}
psSphere *psAberration(psSphere *apparent, const psSphere *actual, const psSphere *direction, double speed);
\end{prototype}
The \code{actual} and \code{apparent} positions are represented as
\code{psSphere} entries, as is the \code{direction} of motion.  The
speed in that direction is given in units of the speed of light.  If
the value of \code{apparent} is NULL, a new \code{psSphere} is
allocated, otherwise the point to \code{apparent} is used for the
result.

\paragraph{Gravitational Deflection}

The following function calculates the \code{apparent} position of a
star, given its \code{actual} position and the position of the sun:
\begin{prototype}
psSphere *psGravityDeflection(psSphere *apparent, psSphere *actual, psSphere *sun);
\end{prototype}
The \code{actual} and \code{apparent} positions are represented as
\code{psSphere} entries, as is position of the sun.  If the value of
\code{apparent} is NULL, a new \code{psSphere} is allocated, otherwise
the point to \code{apparent} is used for the result.

\paragraph{Parallax}

\begin{prototype}
double psEOC_ParallaxFactor(const psSphere *coords, const psTime *time);
\end{prototype}
Calculate the parallax factor for the given position and time.

\subsubsection{Transformation from GCRS to ITRS}

The following functions calculate the components, $X$, $Y$, and $s$,
representing the location of the earth's pole at any moment, or they
determine the velocity of the pole $X'$, $Y'$, $s'$.  We use the
following structure to carry the polar coordinate information.  This
representation may be converted to a rotation between the frames.

\begin{datatype}
typedef struct {
  double x;
  double y;
  double s;
} psEarthPole;
\end{datatype}

\paragraph{Precession/Nutation}

The following routine calculates the components of the rotation
between the CEO and GCRS frames, $X$, $Y$, and $s$, using to the
IAU2000A precession \& nutation model:
%
\begin{prototype}
psEarthPole *psEOC_PrecessionModel(const psTime *time)
\end{prototype}
%
The input to this function is the desired \code{time}, which may be
represented in any format other than UT1.  This routine must give
results identical to the IERS XYS2000A subroutine within the limits of
machine accuracy.

The following function provides interpolated corrections to the $X$
and $Y$ components of the polar coordinates from the tables provided
by the IERS, just as it does for UT1 and polar motion.

\begin{prototype}
psEarthPole *psEOC_PrecessionCorr(const psTime *time, psTimeBulletin bulletin);
\end{prototype}

The polar correction is applied to the $X$ and $Y$ elements of the
rotation to provide higher accuracy.  The spherical rotation term is
generated by providing the polar coordinate to the following function:
\begin{prototype}
psSphereRot *psSphereRot_CEOtoGCRS(const psEarthPole *pole)
\end{prototype}
This function constructs the rotation element as described in the ADD (
The resulting \code{psSphereRot} may be used to determine the rotation
from CIP/CEO to GCRS.  This function must give results identical to
the IERS BPN2000, within the limits of machine accuracy.

\paragraph{Earth Rotation}

The following routine calculates the rotation of the Earth about the CIP:
\begin{prototype}
psSphereRot *psSphereRot_TEOtoCEO(const psTime *time)
\end{prototype}
The IERS code to create the comparable rotation is embedded in
T2C2000, with the Earth Rotation Angle calculated by ERA2000.

\paragraph{Polar Motion}

The following function provides interpolated values of the polar
motion components, $x_p$ and $y_p$, extracted from the IERS tables.  
\begin{prototype}
psEarthPole *psEOC_GetPolarMotion(const psTime *time, psTimeBulletin bulletin);
\end{prototype}

The following function provides tidal corrections to the polar motion
components, $x_p$ and $y_p$, using the Ray model of Simon et al (see
ADD).
\begin{prototype}
psEarthPole *psEOC_PolarTideCorr(const psTime *time);
\end{prototype}

The following function provides the additional corrections due to nutation
terms with periods less than or equal to two days, as well as the
correction to the $s'$ component of the polar motion:
\begin{prototype}
psEarthPole *psEOC_NutationCorr(psTime *time);
\end{prototype}

The following function converts the polar motion corrections to a
spherical rotation using the prescription in the ADD:
\begin{prototype}
psSphereRot *psSphereRot_ITRStoTEO(const psEarthPole *motion);
\end{prototype}
This function should give identical results to the IERS POM2000
subroutine.

\subsubsection{Earth Orientation Wrappers}

The following function generates the complete spherical rotation to
account for precession between two times.  If \code{NULL} is provided
for either time, it is assumed to have the reference equinox value of
J2000.
\begin{prototype}
psSphere *psSpherePrecess(const psTime *fromTime, const psTime *toTime, psPrecessMethod mode);
\end{prototype}
The mode argument is used to specify the level of detail used in the
calculation.

\begin{datatype}
typedef enum {
  PS_PRECESS_ROUGH,
  PS_PRECESS_COMPLETE,
  PS_PRECESS_IAU2000A,
} psPrecessMethod;
\end{datatype}

\subsection{Atmospheric Effects}

\tbd{The ATM effects components should be deferred until we clean up
  the refraction definitions}

A-priori astrometric transformations between the telescope orientation
(Alt/Az) and the predicted stellar coordinates above the atmosphere
(DEC/HA) requires several pieces of information describing the current
environmental conditions.  These quantities are consistent across an
image, and may vary only slowly with time.  Pre-computing these
quantities for exposures means that subsequent transformations are
faster.  The structure below carries the environmental data of interest.
For historical reasons, this structure is known colloquially as ``the
Grommit''.

\tbd{this structure needs to be modified to correspond to what we
  actually need to carry around for the atmosphere functions}

\tbd{provide a single Grommit to carry around all EOC + ATM
  pre-calculated entries and a separate structure for ATM effect?}

\begin{verbatim} %%% Changed to verbatim to remove from api-delta
typedef struct {
    const double latitude;              ///< geodetic latitude (radians)
    const double longitude;             ///< longitude + ... (radians)
    const double height;                ///< height (HM)
    const double abberationMag;         ///< magnitude of diurnal aberration vector
    const double temperature;           ///< ambient temperature (TDK)
    const double pressure;              ///< pressure (PMB)
    const double humidity;              ///< relative humidity (RH)
    const double wavelength;            ///< wavelength (WL)
    const double lapseRate;             ///< lapse rate (TLR)
    const double refractA, refractB;    ///< refraction constants A and B (radians)
    const double siderealTime;          ///< local apparent sidereal time (radians)
} psGrommit;
\end{verbatim}

The \code{psGrommit} is calculated from telescope information for the
particular exposure, \code{exp}:
\begin{verbatim}
psGrommit *psGrommitAlloc(const psExposure *exp);
\end{verbatim}

\tbd{these functions probably need to take the ATM structure}

We require additional functions to perform general functions which
will be useful for astrometry.  Given coordinates on the sky, we
need to get the airmass, the parallactic angle, and an estimate of
the atmospheric refraction.

\begin{verbatim}
float psGetAirmass(const psSphere *coord, psTime *lst, float altitude);
\end{verbatim}
which returns the airmass for a given position and local sidereal time
(\code{lst}).

\begin{verbatim}
float psGetParallactic(const psSphere *coord, double siderealTime);
\end{verbatim}
which returns the parallactic angle for a given position and sidereal time.

\begin{verbatim}
float psGetRefraction(float colour,            ///< Colour of object
                      psPhotSystem colorPlus,  ///< Colour reference
                      psPhotSystem colorMinus, ///< Colour reference
                      const psExposure *exp);  ///< Telescope pointing information
\end{verbatim}
which provides an estimate of the atmospheric refraction, along the parallactic angle.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

\subsection{Projections}

We require functions to convert between spherical and linear
coordinate systems based on a variety of projections.  The required
projections include:
\begin{itemize}
\item TAN
\item SIN
\item AIT
\item PAR
\end{itemize}

We specify the following structure \code{psProjection} to define the
parameters of the projection:
\begin{datatype}
typedef struct {
    double R;                           ///< coordinates of projection center
    double D;                           ///< coordinates of projection center
    double Xs;                          ///< plate-scale in X and Y directions
    double Ys;                          ///< plate-scale in X and Y directions
    psProjectionType type;              ///< projection type
} psProjection;
\end{datatype}

The projection type is defined by the following enumerated type \code{psProjectionType}:
\begin{datatype}
typedef enum {                          ///< type of val is:
    PS_PROJ_TAN,                        ///< Tangent projection
    PS_PROJ_SIN,                        ///< Sine projection
    PS_PROJ_AIT,                        ///< Aitoff projection
    PS_PROJ_PAR,                        ///< Par projection
    PS_PROJ_NTYPE                       ///< Number of types; must be last
} psProjectionType;
\end{datatype}

The constructor is straight-forward:
\begin{prototype}
psProjection *psProjectionAlloc(double R, double D, double Xs, double Ys, psProjectionType type);
\end{prototype}

The following functions will project and deproject (respectively)
spherical coordinates:

\begin{prototype}
psPlane  *psProject(const psSphere *coord, const psProjection *projection);
psSphere *psDeproject(const psPlane *coord, const psProjection *projection);
\end{prototype}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

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

\subsection{Astronomical objects}

\tbd{These functions are all of low priority, have not yet been
defined in detail, and hence are to be deferred.}

\subsubsection{Positions of Major SS Objects}

We may require the ability to calculate the position of major Solar System
objects, as well as Lunar phase.

\begin{verbatim} %%% XXX: This is set to 'verbatim' instead of 'prototype'
psSphere *psSunGetPos(psTime *time);
psTime *psSunGetRise (psTime *twi15, psTime *twi18, const psTime *time);
psTime *psSunGetSet (psTime *twi15, psTime *twi18, const psTime *time);

psSphere *psMoonGetPos(psTime time, psSphere location);
psTime *psMoonGetRise (psTime *twi15, psTime *twi18, psTime *time);
psTime *psMoonGetSet (psTime *twi15, psTime *twi18, psTime *time);
float psGetMoonPhase(psTime time);

psSphere *psPlanetGetPos(psTime time, psSphere location);
psTime *psPlanetGetRise (psTime *twi15, psTime *twi18, psTime *time);
psTime *psPlanetGetSet (psTime *twi15, psTime *twi18, psTime *time);
\end{verbatim}

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

\appendix

\pagebreak
\section{Configuration File Test Inputs}
\label{sec:configtest}

Here are a series of test inputs for the Configuration File syntax defined in
\S\ref{sec:configspec}.

\input{configFileTests.tex}

\section{Dates \& Times Test Inputs}
\label{sec:timetest}

Here are a series of test inputs for the \code{psTime}(\S\ref{sec:timespec})
ISO8601 parsing and formatting functions.  These tests will also validate the
behavior of \code{psTime}'s conversion algorithms.

\input{timeTests.tex}

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

\end{document}
