Index: /tags/ipp-docs/SDRS-08/doc/pslib/psLibSDRS.tex =================================================================== --- /tags/ipp-docs/SDRS-08/doc/pslib/psLibSDRS.tex (revision 22186) +++ /tags/ipp-docs/SDRS-08/doc/pslib/psLibSDRS.tex (revision 22186) @@ -0,0 +1,4749 @@ +%%% $Id: psLibSDRS.tex,v 1.137 2004-10-13 02:34:41 eugene Exp $ +\documentclass[panstarrs,spec]{panstarrs} + +% basic document variables +\title{Pan-STARRS Image Processing Pipeline Library} +\subtitle{Supplementary Design Requirements} +\author{Eugene Magnier, Paul Price, Robert Lupton, Joshua Hoblitt} +\shorttitle{PSLib SDRS} +\group{Pan-STARRS Algorithm Group} +\project{Pan-STARRS Image Processing Pipeline} +\organization{Institute for Astronomy} +\version{08} +\docnumber{PSDC-430-007} +% note the use of the docnumber & version number: +% the complete PSDC document number is given by +% \thedocnumber-\theversion + +\setlength{\topsep}{-2pt} + +\begin{document} +\maketitle +\sloppy + +% -- Revision History -- +% provide explicit values for the old versions +% use '\theversion' for the current version (set above) +% use \hline between each table row +\RevisionsStart +% version Date Description +DR & 2004 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 \\ +\RevisionsEnd + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\DocumentsInternal +PSDC-130-001 & PS-1 Design Reference Mission \\ \hline +PSDC-430-004 & Pan-STARRS IPP C Code Conventions \\ \hline +PSDC-430-005 & Pan-STARRS IPP SRS \\ \hline +PSDC-430-006 & Pan-STARRS IPP ADD \\ \hline +PSDC-430-008 & Pan-STARRS IPP Architecture SDR \\ +\DocumentsExternal +Posix Standard & Open Group Based Specifications Issue 6, IEEE Std 1003.1, 2003 \\ +\DocumentsEnd + +\tableofcontents +\pagebreak +\pagenumbering{arabic} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Introduction} + +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 } 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} + +\item Many of the astronomy routines will wrap the StarLink Positional +Astronomy libraries (SLALib): + +\href{star-www.rl.ac.uk/star/docs/sun67.htx/sun67.html}{\tt +star-www.rl.ac.uk/star/docs/sun67.htx/sun67.html} + +\item libTAI will be used for time-related functions: + +\href{http://cr.yp.to/libtai.html}{\tt http://cr.yp.to/libtai.html} + +\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 the data containers +(doubly-linked lists, hashes, arrays). + +\subsection{Threads} + +Pan-STARRS does not have a strong requirement for multithreading +capability. The memory management functions, defined below, must be +written to be thread-safe. + +\subsection{Conventions} + +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. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\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{verbatim} +const char *psLibVersion(void); +\end{verbatim} + +\subsection{Initialization} + +PSLib shall be initialized by calling \code{psLibInit} by the user +before use. This enables the library to perform the following tasks +that are required before general use: +\begin{itemize} +\item Set the error codes for PSLib. +\item Read the \code{psTime} configuration file (\code{timeConfig}) + and set up the appropriate \code{psTimeTable}s and predictions. +\end{itemize} + +The prototype is: +\begin{verbatim} +void psLibInit(const char *timeConfig); +\end{verbatim} + +\code{timeConfig} specifies the filename of the configuration file +for \code{psTime} metadata. + + +\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 given 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} + +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: +% +\begin{verbatim} +typedef struct { + const void* startblock; ///< initialised to p_psMEMMAGIC + struct psMemBlock* previousBlock; ///< previous block in allocation list + struct psMemBlock* nextBlock; ///< next block allocation list + psFreeFcn freeFcn; ///< deallocator. If NULL, use generic deallocation. + size_t userMemorySize; ///< the size of the user-portion of the memory block + const psMemoryId id; ///< a unique ID for this allocation + const char* file; ///< set from __FILE__ in e.g. p_psAlloc + const 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 (*psFreeFcn)(void* ptr); +typedef unsigned long psMemoryId; +typedef unsigned long psReferenceCount; +\end{verbatim} +% +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{freeFcn} specifies the deallocator associated with a +specific block of memory. If this element is \code{NULL}, the basic +deallocator is used and the memory block must not be a rich data +structure which requires additional freeing functionality. + +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}. + +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. + +\subsubsection{APIs for Allocating and Freeing} + +PSLib must provide the following APIs to create and destroy memory +blocks: +% +\begin{verbatim} +void *psAlloc(size_t size); +void *psRealloc(void *ptr, size_t size); +void *psFree(void *ptr); +\end{verbatim} +% +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{verbatim} +void *p_psAlloc(size_t size, char *file, int line); +void *p_psRealloc(void *ptr, size_t size, char *file, int line); +void p_psFree(void *ptr, char *file, int line); +\end{verbatim} +% +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 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{verbatim} +void p_psMemSetDeallocator(void* ptr, psFreeFcn freeFcn); +psFreeFcn p_psMemGetDeallocator(void* ptr); +\end{verbatim} + +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{p_psSetFreeFcn} is used to associated a free +function to a memory block, while \code{p_psGetFreeFcn} retrieves the +last free function set. To remove a free function from a memory +block, the \code{psSetFreeFcn} should be invoked with \code{NULL}. If +no free function is set, \code{p_psGetFreeFcn} shall return +\code{NULL}. + +\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{psMemAllocateCallbackDefault} +\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{verbatim} +typedef void (*psMemExhaustedCallback)(size_t size); +psMemExhaustedCallback psMemExhaustedCallbackSet(psMemExhaustedCallback func); +\end{verbatim} +% +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{verbatim} +typedef void (*psMemProblemCallback)(psMemBlock *ptr, char *file, int lineno); +psMemProblemCallback psMemProblemCallbackSet(psMemProblemCallback func); +\end{verbatim} +% +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 psMemAllocateCallback \& psMemFreeCallback} + +Two private variables, \code{p_psMemAllocateID} and +\code{p_psMemFreeID}, can be used to trace the allocation and freeing +of specific memory blocks. If the first (\code{p_psMemAllocateID}) 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{verbatim} +psMemoryId psMemAllocateCallbackSetID(psMemoryId id); +psMemoryId psMemFreeCallbackSetID(psMemoryId id); +\end{verbatim} +% +The corresponding callbacks have the following form: +% +\begin{verbatim} +typedef psMemoryId (*psMemAllocateCallback)(const psMemBlock *ptr); +psMemAllocateCallback psMemAllocateCallbackSet(psMemAllocateCallback func); +typedef psMemoryId (*psMemFreeCallback)(const psMemBlock *ptr); +psMemFreeCallback psMemFreeCallbackSet(psMemFreeCallback func); +psMemoryId psMemGetId(void); +\end{verbatim} +% +The callback functions are called with a pointer to the corresponding +memory block. The routines \code{psMemFreeCallbackIDSet} and +\code{psMemAllocateCallbackIDSet} accept the desired ID value and return the +old value to the user. The return values of the handlers installed by +\code{psMemAllocateCallbackSet} and \code{psMemFreeCallbackSet} are +used to increment the values of \code{p_psMemAllocateID} 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{verbatim} +int psMemCheckLeaks(psMemoryId id0, psMemBlock ***array, FILE *fd, bool persistence); +int psMemCheckCorruption(bool abort_on_error); +\end{verbatim} +% +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{verbatim} +psReferenceCount psMemGetRefCounter(void *vptr); +void *psMemIncrRefCounter(void *vptr); +void *psMemDecrRefCounter(void *vptr); +\end{verbatim} +% +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{psMyType} should be constructed using +\code{psMyTypeAlloc()} calls, and are destroyed using the basic +\code{psFree} function, which calls \code{psMyTypeFree()} 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 +(i.e., \code{psMyTypeFree}) and reference counters. Simply put, +\textit{the destructor for every structure should only call +\code{psFree} if \code{refCounter == 1}; otherwise, it decrements the +reference counter and returns.} An example destructor is shown below: + +\begin{verbatim} +void psMyTypeFree(psMyType *myType) +{ + /* data is not defined */ + if (psMemGetRefCounter(myType) < 1) { + return; + } + /* Only call psFree if reference counter is 1 */ + if (psMemGetRefCounter(myType) == 1) { + psSubFree (myType->sub); + psFree(myType); + return; + } + /* Otherwise, decrement the reference counter only */ + psMemDecrRefCounter(myType); +} +\end{verbatim} + +Note that the element of \code{myType}, \code{myType.sub} is +explicitly freed with its associated destructor. If this element +points to a data block with multiple references, this call would only +decrement the counter. + +\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{verbatim} +void psTrace(char *facil, int myLevel,...); +\end{verbatim} +% +where the last argument is a printf-style formatting code and 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{myLevel}) and the message +itself. + +The trace level for any facility may be set at any time with the +following function: +% +\begin{verbatim} +int psTraceSetLevel(const char *facil, int level); +\end{verbatim} +% +where \code{level} specifies the current trace level for the facility +named by \code{facil}. The currently defined trace level for a given +facility may be determined by the function: +% +\begin{verbatim} +int psTraceGetLevel(const char *facil); +\end{verbatim} +% +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{myLevel} 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{verbatim} +void psTracePrintLevels(void); +\end{verbatim}. +% +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{void psTraceReset(void)} will set all trace levels to 0. + + +The trace may optionally be written to a file on specifying a filename +with \code{psTraceSetDestination}: +\begin{verbatim} +void psTraceSetDestination(FILE *fp); +\end{verbatim} +If the \code{fp} is \code{NULL}, then the trace is sent to standard +output, otherwise it is sent to the specified file pointer. + +\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{verbatim} +void psLogMsg(const char *name, int myLevel, const char *fmt, ...); +void psLogMsgV(const char *name, int myLevel, const char *fmt, va_list ap); +\end{verbatim} +where \code{name} is a word to describe the source of the message, +\code{myLevel} is the severity level of this message, and \code{fmt} +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{verbatim} +enum { PS_LOG_ABORT = 0, PS_LOG_ERROR, PS_LOG_WARN, PS_LOG_INFO }; +\end{verbatim} +% + +At any time, the program may set the current log level, the level +above which log messages are ignored, using the function: +% +\begin{verbatim} +int psLogSetLevel(int level); +\end{verbatim} +% +This function returns the previous log level. A specific message +invoked with \code{psLogMsg} is only printed if its value of +\code{myLevel} is less than the current value set by +\code{psLogSetLevel}. The default log level is set to +\code{PS_LOG_INFO}. + +Log messages are sent to the destination most recently set using: +% +\begin{verbatim} +bool psLogSetDestination(const char *dest); +\end{verbatim} +% +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 write to \code{stderr} and \code{stdout} respectively, and +\code{none} to turn off logging. + +The output format is controlled with the function: +% +\begin{verbatim} +void psLogSetFormat(const char *fmt); +\end{verbatim} +% +which expects a string consisting of the letters \code{H} (host), +\code{L} (level), \code{M} (message), \code{N} (name), 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 + 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{fmt} of \code{NULL} turns off logging (equivalent +to calling \code{psLogSetDestination(PS_LOG_TO_NONE)}, whereas if the +\code{fmt} is \code{""}, then the format reverts to the default. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Error Handling} +\hlabel{errorStack} + +\PS{} errors must 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{verbatim} +psErrorCode p_psError(const char *file, int lineno, psErrorCode code, bool new, const char *fmt, ...); +#define psError(code, new, fmt, ...) psError(__FILE__, __LINE__, code, new, fmt, __VA_ARGS__) +\end{verbatim} + +\code{psError} is a macro definition that allows the function and line +number 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{fmt}, 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} must 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{verbatim} +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{verbatim} + +The last error reported is available from \code{psErrorLast}; if no +errors are current, a non-\code{NULL} \code{psErr} must be returned +with code \code{PS_ERR_NONE}. Previous errors on the stack must 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{verbatim} +const psErr *psErrorGet(int which); +const psErr *psErrorLast(void); +void psErrorClear(void); +\end{verbatim} + +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{fmt} is first printed to the file descriptor \code{fd}. In this +printout, error codes must 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{verbatim} +void psErrorStackPrint(FILE *fd, const char *fmt, ...); +void psErrorStackPrintV(FILE *fd, const char *fmt, va_list va); +\end{verbatim} + +Any \code{errorCode}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{verbatim} +const char *psErrorCodeString(psErrorCode code); +\end{verbatim} + +\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, 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{verbatim} +void psErrorRegister(const psErrorDescription *errors, int nerror); +\end{verbatim} +where the errors are represented internally as follows: +\begin{verbatim} +typedef struct { + psErrorCode code; ///< An error code + const char *descrip; ///< the associated description +} psErrorDescription; +\end{verbatim} +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{verbatim} +#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{verbatim} + +Any \code{errorCode}s less then or equal to \code{PS_ERR_BASE} must be +taken to be valid values of \code{errno}. \file{psErrorCodes.c} must +define the necessary function to register the error codes. + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Abort} + +\code{psAbort}, must call \code{psLogMsgV} with a level of +\code{PS_LOG_ABORT}, and then call \code{abort}. + +\begin{verbatim} +void psAbort(char *name, char *fmt,...); +\end{verbatim} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Basic Data Types and Collections} + +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; and +\item Hashes. +\end{itemize} + +\subsection{Data Structure Type Information} + +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{psType}: +\begin{verbatim} +typedef struct { + psDimen dimen; ///< The dimensionality + psElemType type; ///< The type +} psType; +\end{verbatim} +where \code{psDimen dimen} defines the dimensionality of the data and +\code{psElemType type} defines the data type of each element. These +two variable types are defined as: +\begin{verbatim} +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{verbatim} +and +\begin{verbatim} +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{verbatim} +We discuss the application of \code{psType} in more detail in +section~\ref{sec:arithmetic}. + +\subsection{Simple 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{verbatim} +typedef struct { + psType type; ///< data type information + union { + psS8 S8; ///< bye value entry + psS16 S16; ///< short int value entry + psU8 U8; ///< unsigned byte value entry + psU16 U16; ///< unsigned short 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{verbatim} + +In addition, we specify two functions for working with \code{psScalar} data: +\begin{verbatim} +psScalar *psScalarAlloc (psC64 value, psElemType dataType); +psScalar *psScalarCopy (psScalar *value); +\end{verbatim} +The first creates a \code{psType}-ed structure from a constant value, +casting it as appropriate based on the \code{dataType}. 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{Simple 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 double}. We have defined a +single structure, \code{psVector} to represent these concepts: +% +\begin{verbatim} +typedef struct { + psType type; ///< vector data type and dimension + const int n; ///< size of vector + const int 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 *U18; ///< 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; +} psVector; +\end{verbatim} +% +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 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 psElemType names. This parallelism allows us +to use automatic construction mechanisms effectively. The data type +is defined by the first element, \code{psType}. The structure is +associated with a constructor and reallocator: +% +\begin{verbatim} +psVector *psVectorAlloc(int nalloc, psElemType type); +psVector *psVectorRealloc(const psVector *vector, int nalloc); +\end{verbatim} +% +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 return an error. + +\subsection{Simple 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{verbatim} +typedef struct psImage { + psType type; ///< image data type and dimension + const int ncols, nrows; ///< size of image + const int x0, y0; ///< data region 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 + } data; + const struct psImage *parent; ///< parent, if a subimage + psArray *children; ///< children of this region +} psImage; +\end{verbatim} + +This structure represents an image consisting of a 2-D array of +pixels. The size of this array is given by the elements \code{(nrows, +ncols)}. The data type of the pixel is defined in the \code{psType +type} entry (specifically, the \code{psElemType} member, \code{type}; +see \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{(x0,y0)}: \code{x0} is the +starting column number in the parent image while \code{y0} is the +starting row number. The structure may include references to +subrasters (\code{children, nChildren}) 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. We define the following +supporting functions, which are valid for data types \code{psS8, +psS16, psU8, psU16, psF32, psF64, psC32, psC64}. + +\begin{verbatim} +psImage *psImageAlloc (int width, int height, psElemType type); +\end{verbatim} +Create an image of a specified \code{width}, \code{height}, 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. + +\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{verbatim} +typedef struct { + psType type; ///< data type information +} psMath; +\end{verbatim} + +\subsection{Simple Arrays} + +We require an order collection of unspecified data elements. We +define \code{psArray} to carry such a collection: +% +\begin{verbatim} +typedef struct { + const int n; ///< size of array + const int nalloc; ///< allocated data block + void **data; ///< pointer to data block +} psArray; +\end{verbatim} +% +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{verbatim} +psArray *psArrayAlloc(int nalloc); +psArray *psArrayRealloc(const psArray *array, int nalloc); +\end{verbatim} +% +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{verbatim} +psArray *psArraySort(psArray *array, int (*compare)(const void **a, const void **b) ); +\end{verbatim} +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. + +\subsection{Doubly-linked lists} +\label{sec:psList} + +\PS{} shall support doubly linked lists through a type \code{psList}: +% +\begin{verbatim} +typedef struct { + unsigned int size; ///< number of elements on list + psListElem *head; ///< first element on list (may be NULL) + psListElem *tail; ///< last element on list (may be NULL) + unsigned int nIter; ///< number of iterators + psListIterator **iter; ///< iteration cursors + pthread_mutex_t lock; ///< mutex to lock a node during changes +} psList; +\end{verbatim} +% +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 entry for +the current cursor location (\code{iter}), and an entry to define the +number of elements in the list (\code{size}). Iteration on the list +is achieved by means of an array of \code{nIter} iteration cursors, +\code{iter}. The \code{lock} should be applied during changes to the +list to maintain thread safety. + +The elements of the list are defined by the type \code{psListElem}: +% +\begin{verbatim} +typedef struct psListElem { + struct psListElem *prev; ///< previous link in list + struct psListElem *next; ///< next link in list + void *data; ///< real data item +} psListElem; +\end{verbatim} +% +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{verbatim} +psList *psListAlloc(void *data); +\end{verbatim} +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 are set to +\code{NULL} should be created. The number of iterators in the list is +initially set to zero. + +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) also increments the reference +counter. + +Iteration on the list shall be achieved by means of a list iterator +type: +\begin{verbatim} +typedef struct { + psList *list; ///< List iterator works on + unsigned int num; ///< Iterator number + psListElem *cursor; ///< The current iterator cursor + bool offEnd; ///< Is the iterator off the end? +} psListIterator; +\end{verbatim} +The \code{psListIterator} keeps track of which list element the +iterator \code{cursor} is currently pointing at. \code{num} is the +number of the list iterator in the parent \code{list}'s array of +iterators, so \code{list->iter[num]} points to the list iterator. The +boolean member, \code{offEnd}, indicates whether the iterator has +progressed off the end of the list (i.e., beyond the last item). + +The corresponding constructor shall be: +\begin{verbatim} +psListIterator *psListIteratorAlloc(psList *list, int location); +\end{verbatim} +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 +negative integers. Of course, the \code{psListElem} pointed to by the +iterator shall have its \code{refCounter} incremented. + +The destructor for \code{psListIterator} shall, after freeing the +\code{psListIterator} and decrementing the \code{refCounter} of the +\code{psListElem} that it pointed to, 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{verbatim} +bool psListIteratorSet(psListIterator *iterator, int location); +\end{verbatim} +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 negative integers. The function shall return +\code{true} if the reset was successful, or \code{false} otherwise. + +\begin{verbatim} +bool psListAdd(psList *list, int location, void *data); +bool psListAddAfter(psListIterator *iterator, void *data); +bool psListAddBefore(psListIterator *iterator, void *data); +\end{verbatim} +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 negative integers. 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{verbatim} +void *psListGet(psList *list, int location); +void *psListGetNext(psListIterator *iterator); +void *psListGetPrevious(psListIterator *iterator); +\end{verbatim} +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} or \code{PS_LIST_TAIL}, +which are defined as negative integers, allowing \code{location} to +also be the index of one of the data items. The other two functions, +\code{psListGetNext} and \code{psListGetPrevious} 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{psListGetNext} 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{psListGetPrevious} 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{psListGetNext} and \code{psListGetPrevious} 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{psListGetPrevious} 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. +When setting the \code{cursor}, the appropriate \code{refCounter}s +shall be incremented and decremented. + +\begin{verbatim} +bool psListRemove(psList *list, int location) +bool psListRemoveData(psList *list, void *data); +\end{verbatim} +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 negative integers. For +\code{psListRemoveData}, the data item to be removed is identified by +matching the pointer value with \code{void *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{verbatim} +psArray *psListToArray(psList *list); +psList *psArrayToList(psArray *array); +\end{verbatim} +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{verbatim} +psList *psListSort(psList *list, int (*compare)(const void **a, const void **b) ); +\end{verbatim} +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. + +\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{verbatim} +typedef struct { + int nbucket; ///< number of buckets + psHashBucket **buckets; ///< the buckets themselves +} psHash; +\end{verbatim} +% +where \code{nbucket} 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{verbatim} +typedef struct psHashBucket { + char *key; ///< key for this item of data + void *data; ///< the data itself + struct psHashBucket *next; ///< list of other possible keys +} psHashBucket; +\end{verbatim} +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{verbatim} +psHash *psHashAlloc(void); +\end{verbatim} +which allocates the space for the hash table and initializes all of +the 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{verbatim} +bool psHashAdd(psHash *table, const char *key, void *data); +\end{verbatim} +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{verbatim} +void *psHashLookup(psHash *table, const char *key); +\end{verbatim} +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{verbatim} +bool psHashRemove(psHash *table, const char *key); +\end{verbatim} +The function returns a value of \code{true} if the operation was +successfull, and \code{false} otherwise. + +The function +\begin{verbatim} +psList *psHashKeyList(psHash *table); +\end{verbatim} +returns the complete list of defined keys associated with the +\code{psHash} table as a linked list. + +\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 Bit masks; +\item Vector and image arithmetic; +\item Sorting; +\item Statistics; +\item Matrix operations and linear algebra; +\item (Fast) Fourier Transforms; +\item General functions; and +\item Minimization and fitting routines. +\end{itemize} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\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{verbatim} +typedef struct { + int n; ///< Number of chars that form the bitset + char *bits; ///< The bits +} psBitSet; +\end{verbatim} + +We also require the corresponding constructor and destructor: +\begin{verbatim} +psBitSet *psBitSetAlloc(int n); +\end{verbatim} +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{verbatim} +psBitSet *psBitSetSet(psBitSet *myBits, int bit); +psBitSet *psBitSetOp(psBitSet *outBits, + const psBitSet *inBits1, + char *operator, + const psBitSet *inBits2); +psBitSet *psBitSetNot(psBitSet *out, psBitSet *in); +bool psBitSetTest(const psBitSet *checkBits, int bit); +char *psBitSetToString(const pBitSet* bits); +\end{verbatim} + +\code{psBitSetSet} sets the specified \code{bit} in the +\code{psBitSet}, 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{inBits1} and \code{inBits2}. +If the output bitset \code{outBits} 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{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{verbatim} +psVector *psVectorSort(psVector *out, const psVector *in); +\end{verbatim} + +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{verbatim} +psVector *psVectorSortIndex(psVector *out; const psVector *in); +\end{verbatim} + +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{Statistics Functions} + +\subsubsection{Vector Statistics} + +We require a very general statistics function, which, given a vector +of floating-point values, will be able to calculate the following +statistics: +\begin{itemize} +\item Sample mean; +\item Sample median; +\item Sample standard deviation; +\item Sample upper and lower quartiles; +\item Robust mean and number of values used to calculate; +\item Robust median and number of values used to calculate; +\item Robust mode and number of values used to calculate; +\item Robust standard deviation; +\item Robust upper and lower quartiles; +\item Clipped mean and number of values used to calculate; +\item Clipped standard deviation; and +\item Minimum and maximum value in vector. +\end{itemize} + +For definitions of each of these, see the accompanying Algorithms +Definition Document (ADD), but in general, ``sample'' refers to the +entire vector, ``robust'' refers to fitting the distribution in the +vector with a model, and ``clipped'' refers to clipping the +distribution. Each of these must be available from a single +function: + +\begin{verbatim} +psStats *psVectorStats(psStats *stats, + const psVector *in, + const psVector *errors, + const psVector *mask, + unsigned int maskVal + ); +\end{verbatim} +% +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{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 mask must be of type \code{psU8}. + +The \code{psStats} structure is defined with entries for each of the +desired statistical quantities: + +\begin{verbatim} +/** generic statistics structure */ +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 robustMean; ///< robust mean of data + double robustMedian; ///< robust median of data + double robustMode; ///< Robust mode 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 + int robustNfit; ///< 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) + psStatsOptions options; ///< bitmask of calculated values +} psStats; +\end{verbatim} +where \code{psStatsOptions} is defined with entries to turn on the +calculation of each of the statistics: + +\begin{verbatim} +/** statistics which may be calculated */ +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_MEAN = 0x000010, + PS_STAT_ROBUST_MEDIAN = 0x000020, + PS_STAT_ROBUST_MODE = 0x000040, + PS_STAT_ROBUST_STDEV = 0x000080, + PS_STAT_ROBUST_QUARTILE = 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{verbatim} + +A constructor is also required: +% +\begin{verbatim} +psStats *psStatsAlloc(psStatsOptions options); +\end{verbatim} + +\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{verbatim} +typedef struct { + const psVector *bounds; ///< Bounds for the bins + psVector *nums; ///< Number in each of the bins + int minNum, maxNum; ///< Number below minimum / above maximum + bool uniform; ///< Is it a uniform distribution? +} psHistogram; +\end{verbatim} +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 of +type \code{psU32}. 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{verbatim} +psHistogram *psHistogramAlloc(float lower, float upper, int n); +\end{verbatim} +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{verbatim} +psHistogram *psHistogramAllocGeneric(const psVector *bounds); +\end{verbatim} +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{in}), and optionally the \code{errors} in the input +values. It alters and returns the histogram \code{out} structure. +\begin{verbatim} +psHistogram *psVectorHistogram(psHistogram *out, + const psVector *in, + const psVector *mask, + unsigned int maskVal); +\end{verbatim} +The input vector may be of types \code{psU8, psU16, psF32, psF64}. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\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. + +This leads us to define the following polynomial types: + +\begin{verbatim} +/** One-dimensional polynomial */ +typedef struct { + psPolynomialType type; ///< Polynomial type + int n; ///< Number of terms + float *restrict coeff; ///< Coefficients + float *restrict coeffErr; ///< Error in coefficients + char *restrict mask; ///< Coefficient mask +} psPolynomial1D; +\end{verbatim} + +\begin{verbatim} +/** Two-dimensional polynomial */ +typedef struct { + psPolynomialType type; ///< Polynomial type + int nX, nY; ///< Number of terms in x and y + float *restrict *restrict coeff; ///< Coefficients + float *restrict *restrict coeffErr; ///< Error in coefficients + char *restrict *restrict mask; ///< Coefficients mask +} psPolynomial2D; +\end{verbatim} + +etc., up to four dimensions. \code{psPolynomialType} is an +enumerated type specifying the type of the polynomial: ordinary +or Chebyshev: + +\begin{verbatim} +typedef enum { + PS_POLYNOMIAL_ORD, ///< Ordinary polynomial + PS_POLYNOMIAL_CHEB ///< Chebyshev polynomial +} psPolynomialType; +\end{verbatim} + +We also define double-precision versions of the polynomials: + +\begin{verbatim} +/** Double-precision one-dimensional polynomial */ +typedef struct { + psPolynomialType type; ///< Polynomial type + int n; ///< Number of terms + double *restrict coeff; ///< Coefficients + double *restrict coeffErr; ///< Error in coefficients + char *restrict mask; ///< Coefficient mask +} psDPolynomial1D; +\end{verbatim} + +\begin{verbatim} +/** Double-precision two-dimensional polynomial */ +typedef struct { + psPolynomialType type; ///< Polynomial type + int nX, nY; ///< Number of terms in x and y + double *restrict *restrict coeff; ///< Coefficients + double *restrict *restrict coeffErr; ///< Error in coefficients + char *restrict *restrict mask; ///< Coefficients mask +} psDPolynomial2D; +\end{verbatim} + +etc. In what follows, we only show the version for double-precision +two-dimensionals; the others may be inferred following the standard +naming convention exampled above. + +The constructor is: +\begin{verbatim} +psDPolynomial2D *psDPolynomial2DAlloc(int nX, int nY, psPolynomialType type); +\end{verbatim} +where \code{nX} and \code{nY} are the number of terms in x and y +respectively. The coefficients and errors are set initially to 0.0. + +To evaluate the polynomials at specific coordinates, we define: +\begin{verbatim} +double psDPolynomial2DEval(double x, double y, const psDPolynomial2D *restrict myPoly); +\end{verbatim} + +In the event that several evaluations are required, we also define: +\begin{verbatim} +psVector *psDPolynomial2DEvalVector(const psVector *x, const psVector *y, + const psDPolynomial2D *myPoly); +\end{verbatim} +If the \code{x} and \code{y} vectors do not match the precision of the +polynomial, the function shall generate a warning, and convert the +vectors to the appropriate precision. 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 continue. The precision of the output \code{psVector} +shall match that of the polynomial. + +\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{verbatim} +typedef struct { + int n; ///< Number of spline pieces + psPolynomial1D **spline; ///< Array of n pointers to splines + float *domains; ///< The boundaries between each spline piece. Size is n+1. +} psSpline1D; +\end{verbatim} + +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{domains} member specifies the boundaries +between each spline piece (including the two ends). + +Of course, we require the appropriate constructors and destructor: +\begin{verbatim} +psSpline1D *psSpline1DAlloc(int n, int order, float min, float max); +psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, int order); +\end{verbatim} + +\code{psSpline1DAlloc} shall allocate and return a \code{psSpline1D}, +setting the \code{domains} 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{domains}. 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{verbatim} +float psSpline1DEval(const psSpline1D *spline, float x); +psVector *psSpline1DEvalVector(const psVector *x, const psSpline1D *spline); +\end{verbatim} + + +\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{verbatim} +float psGaussian(float x, float mean, float sigma, bool normal); +\end{verbatim} +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 always obtain the best fit but does +not require derivatives to be known. + +In both cases, we need to be able to set the maximum number of +iterations (\code{maxIter}) and the maximum tolerance (\code{tol}) for +convergence, in addition to having some basic information on the +results and quality of the minimization: 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{verbatim} +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 +} psMinimization; +\end{verbatim} + +The corresponding allocator is: +\begin{verbatim} +psMinimization *psMinimizationAlloc(int maxIter, float tol); +\end{verbatim} + +\subsubsection{Levenberg-Marquardt} + +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}, 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: +\begin{verbatim} +typedef float (*psMinimizeLMFunc)(psVector *deriv, const psVector *params, const psArray *coords); +\end{verbatim} +An example may help here. We say that the coordinate vectors may be +several (and hence use a \code{psArray} for \code{coords}), since +often one will want to minimize $\chi^2$, given a data set consisting +of vectors for $x$, $y$ and $\sigma_y$. + +Then \code{psMinimizeLM} shall mimimize the specified function, +\code{func}, using the Levenberg-Marquardt method: + +\begin{verbatim} +bool psMinimizeLM(psMinimization *min, psImage *covar, psVector *params, const psVector *paramMask, + const psArray *coords, psMinimizeLMFunc func); +\end{verbatim} + +The function shall return \code{false} in the event that there was an +error (bad input, too many iterations), and also call \code{psError}. + +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, \code{params} shall consist of a ``best +guess'' for the parameters that minimize the input 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 not a square +matrix with size matching that of \code{params}. + +Parameters that have a corresponding \code{paramMask} entry of +\code{0} 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 \code{coords} are not handled by the minimizer, except to be +passed to the function as additional input, and it is the +responsibility of the function to interpret these. + +\code{paramMask} must be of type \code{psU8}, while \code{params} must +be of type \code{psF32}. The \code{func} function must be valid only +for types \code{psF32}, \code{psF64}. + +\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{verbatim} +psMinimizeLMFunc psMinimizeLMChi2Gauss1D; +psMinimizeLMFunc psMinimizeLMChi2Gauss2D; +\end{verbatim} + +\code{psMinimizeChi2LMGauss1D} shall take, as \code{params}, the +normalization, center, and standard deviation of a Gaussian to be fit +to the data. The data in the \code{coords} shall consist of vectors +for each of \code{x}, \code{y}, and \code{yErr}. If the \code{yErr} +is missing, or \code{NULL}, then the errors shall be taken to be +unity, for the purpose of fitting where the errors are all identical +(and possibly unknown). If the \code{x}, \code{y} or \code{yErr} (if +supplied) vectors are of differing lengths, the function shall +generate a warning and truncate the lengths to match the shortest. + +\code{psMinimizeLMChi2Gauss1D} shall return the $\chi^2$ value for the +specified Gaussian as a model for the supplied data, along with the +derivatives of the $\chi^2$ 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 to be fit to the data. +The data in the \code{coords} shall consist of vectors for each of +\code{x}, \code{y}, \code{z} and \code{zErr}. If the \code{zErr} is +missing, or \code{NULL}, then the errors shall be taken to be unity, +for the purpose of fitting where the errors are all identical (and +possibly unknown). If the \code{x}, \code{y}, \code{z} or \code{zErr} +(if supplied) vectors are of differing lengths, the function shall +generate a warning and truncate the lengths to match the shortest. + +\code{psMinimizeLMChi2Gauss2D} shall return the $\chi^2$ value for the +specified Gaussian as a model for the supplied data, along with the +derivatives of the $\chi^2$ 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{verbatim} +typedef float (*psMinimizePowellFunc)(const psVector *params, const psArray *coords); +\end{verbatim} + +Then \code{psMinimizePowell} shall mimimize the specified function, +\code{func}, using the Powell method: + +\begin{verbatim} +bool psMinimizePowell(psMinimization *min, psVector *params, const psVector *paramMask, + const psArray *coords, psMinimizePowellFunc func); +\end{verbatim} + +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}. + +\subsubsubsection{Minimizing $\chi^2$ with Powell} + +We require a front-end to the Powell minimizer to allow fitting +general functions to data. + +\begin{verbatim} +typedef psVector* (*psMinimizeChi2PowellFunc)(const psVector *params, const psArray *coords); + +bool psMinimizeChi2Powell(psMinimization *min, psVector *params, const psVector *paramMask, + const psArray *coords, const psVector *value, const psVector *error, + psMinimizeChi2PowellFunc model); +\end{verbatim} + +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. + +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{verbatim} +psPolynomial1D *psVectorFitPolynomial1D(psPolynomial1D *myPoly, + const psVector *x, + const psVector *y, + const psVector *yErr); +\end{verbatim} +\code{psVectorFitPolynomial1d} returns the polynomial that best fits the +observations. The input parameters are a polynomial that specifies +the fit order, \code{myPoly}, which will be altered and returned with +the best-fit coefficients; and the observations, \code{x}, \code{y} +and \code{yErr}. 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 types \code{psF32}, \code{psF64}. + +\begin{verbatim} +psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline, + const psVector *x, + const psVector *y, + const psVector *yErr); +\end{verbatim} +\code{psVectorFitSpline} shall return the spline that best fits the +given combination of ordinates (\code{x}) and coordinates (\code{y}). +As for \code{psVectorFitPolynomial1D}, if \code{x} is \code{NULL}, +then the vector index shall be used as the ordinate; and if +\code{yErr} is \code{NULL}, then all the data errors shall be assumed +to be equal. 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{verbatim} +psImage *psImageSubset(psImage *image, int x0, int y0, int x1, int y1); +\end{verbatim} +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{x0, y0} specify the starting pixel of the +subraster, and \code{x1,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{verbatim} +psImage *psImageSubsection(psImage *image, const char *section); +\end{verbatim} +Similar to \code{psImageSubset}, but uses an image \code{section}, of +the form \code{[x1:x2,y1:y2]}. If the image section is not fully +contained within the input \code{image}, then the function shall +generate an error and return \code{NULL}. 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{verbatim} +psImage *psImageCopy(psImage *output, const psImage *input, psElemType type); +\end{verbatim} +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. 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{verbatim} +psImage *psImageTrim(psImage *image, int x0, int y0, int x1, int y1); +\end{verbatim} +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{x0,y0}, and the upper-right corner, +\code{x1,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{x1} or \code{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{verbatim} +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; + +psVector *psImageSlice(psVector *out, psVector *coords, const psImage *input, + const psImage *mask, unsigned int maskVal, int x0, int y0, + int x1, int y1, psImageCutDirection direction, const psStats *stats); +\end{verbatim} +Extract pixels from rectlinear region to a vector (array of floats). +The output vector contains either \code{x1-x0} or \code{y1-y0} +elements, based on the value of the direction: e.g., if +\code{direction} is \code{PS_CUT_X_POS}, there are \code{x1-x0} +elements. The region to be ``sliced'' is defined by the lower-left +corner, \code{x0,y0}, and the upper-right corner, \code{x1,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{x1} or +\code{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 +vector. 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{verbatim} +psVector *psImageCut(psVector *out, + psVector *coords, + const psImage *input, + const psImage *mask, + unsigned int maskVal, + float x0, + float y0, + float x1, + float y1, + unsigned int nSamples, + psImageInterpolateMode mode); +\end{verbatim} +Extract pixels along a line segment, \code{(x0,y0)} to \code{(x1,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{verbatim} +psVector *psImageRadialCut(psVector *out, + const psImage *input, + const psImage *mask, + unsigned int maskVal, + float x, + float y, + const psVector *radii, + const psStats *stats); +\end{verbatim} +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{verbatim} +typedef enum { + PS_INTERPOLATE_FLAT, + PS_INTERPOLATE_BILINEAR, + PS_INTERPOLATE_BICUBIC, + PS_INTERPOLATE_SINC +} psImageInterpolateMode mode; +\end{verbatim} + +\begin{verbatim} +psImage *psImageRebin(psImage *out, const psImage *in, + const psImage *restrict mask, + unsigned int maskVal, + int scale, const psStats *stats); +\end{verbatim} +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{psU8}, +\code{psU16}, \code{psS8}, \code{psS16}, \code{psF32}, \code{psF64}, +\code{psC32}, \code{psC64}. + +\begin{verbatim} +psImage *psImageResample(psImage *out, const psImage *in, + int scale, psImageInterpolateMode mode); +\end{verbatim} +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{verbatim} +psImage *psImageRotate(psImage *out, const psImage *input, float angle, + psC64 exposed, psImageInterpolateMode mode); +\end{verbatim} +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{verbatim} +psImage *psImageShift(psImage *out, const psImage *input, + float dx, float dy, psC64 exposed, psImageInterpolateMode mode); +\end{verbatim} +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{verbatim} +psImage *psImageRoll(psImage *out, const psImage *input, int dx, int dy); +\end{verbatim} +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}. + +\subsubsection{Image Statistical Functions} + +\begin{verbatim} +psStats *psImageStats(psStats *stats, + const psImage *in, + const psImage *mask, + unsigned int maskVal); +\end{verbatim} +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{verbatim} +psHistogram *psImageHistogram(psHistogram *out, + const psImage *in, + const psImage *mask, + unsigned int maskVal); +\end{verbatim} +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{verbatim} +psPolynomial2D *psImageFitPolynomial(psPolynomial2D *coeffs, const psImage *input); +\end{verbatim} +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{verbatim} +psImage *psImageEvalPolynomial(psImage *input, const psPolynomial2D *coeffs); +\end{verbatim} +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{verbatim} +psC64 psImagePixelInterpolate(const psImage *input, float x, float y, + const psImage *mask, unsigned int maskVal, + psC64 unexposedValue, psImageInterpolateMode mode); +\end{verbatim} +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 I/O Functions} + +\tbd{this section is being deprecated and the equivalent of these +functions are being moved to FITS I/O Functions. I will leave this +section here until others may read the section on FITS I/O.} + +\begin{verbatim} +psImage *psImageReadSection(psImage *output, int x0, int y0, int x1, int y1, int z, + const char *extname, int extnum, const char *filename); +\end{verbatim} +Read an image or subimage from a named file. This function is a +wrapper to the FITS library function. The input parameters allow a +full image or a subimage to be read. The starting pixel of the region +is specified by \code{x0,y0}, while the upper right-hand corner of the +requested region is specified by \code{x1,y1}. A negative value for +either of \code{x1} or \code{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. The data is read from the extension specified by +extname (matching the EXTNAME keyword) or by the extnum value (with 0 +representing the primary header unit (PHU), 1 the first extension, +etc). 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, if the specified image file does not exist, or 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{verbatim} +bool psImageWriteSection(const psImage *input, int x, int y, int z, + const char *extname, int extnum, const char *filename); +\end{verbatim} +Write an image section to the named file, which may exist. This +operation may write a portion of an image over the existing bytes of +an existing image. If the file does not exist, it should be created. +If the specified extension does not exist, it should be created. If +an extension is specified and no PHU exists, a basic PHU should be +created. Care must be taken to interpret x,y,z in the two cases a) +there is already an existing image and b) there is not an existing +image. If the image exists, write the complete psImage data to the +existing image starting at the coordinate x,y,z. If any of these +parameters implies writing pixels outside the existing data area of +the image, return an error (ie, if \code{x + image.nx >= NAXIS1}, +\code{y + image.ny >= NAXIS2}, or \code{z >= NAXIS3}). If the image +does not exist, require x,y,z to be zero. 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. + +\subsubsection{Image Pixel Manipulations} + +\begin{verbatim} +int psImageClip(psImage *input, double min, double vmin, double max, double vmax); +\end{verbatim} +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{verbatim} +int psImageClipComplexRegion(psImage *input, complex double min, complex double vmin, + complex double max, complex double vmax); +\end{verbatim} +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{verbatim} +int psImageClipNaN(psImage *input, float value); +\end{verbatim} +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{verbatim} +int psImageOverlaySection(psImage *image, const psImage *overlay, + int x0, int y0, const char *op); +\end{verbatim} +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. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\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{verbatim} +psType *psBinaryOp (void *out, void *in1, char *op, void *in2); +psType *psUnaryOp (void *out, void *in, char *op); +\end{verbatim} +These functions determine the type of the operands on the basis of +their \code{psType} 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{psType}, 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) +\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. + +\begin{verbatim} +psImage *psMatrixLUD(psImage *out, psVector *perm, const psImage *in); +psVector *psMatrixLUSolve(psVector *out, const psImage *LU, const psVector *RHS, const psVector *perm); +\end{verbatim} +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. + +\begin{verbatim} +psImage *psMatrixInvert(psImage *out, const psImage *in, float *determinant); +\end{verbatim} +\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{verbatim} +float psMatrixDeterminant(const psImage *in); +\end{verbatim} + +Matrix multiplication is supported through \code{psMatrixMultiply}. This function is +specified for data types \code{psF32, psF64}. +\begin{verbatim} +psImage *psMatrixMultiply(psImage *out, const psImage *in1, const psImage *in2); +\end{verbatim} + +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{verbatim} +psImage *psMatrixTranspose(psImage *out, const psImage *in); +\end{verbatim} + +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{verbatim} +psVector *psMatrixEigenvectors(psImage *in); +\end{verbatim} + +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, utilised, 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{verbatim} +psVector *psMatrixToVector(psVector *out, const psImage *in); +psImage *psVectorToMatrix(psImage *out, const psVector *in); +\end{verbatim} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\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{verbatim} +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, psVector *real, psVector *imag); +psVector *psVectorConjugate(psVector *out, const psVector *in); +psVector *psVectorPowerSpectrum(psVector *out, const psVector *in); +\end{verbatim} + +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{psVectorConjuagte} 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{verbatim} +/** Specify direction of FFT, and if the result is real or not */ +typedef enum { + PS_FFT_FORWARD = 0, ///< psImageFFT/psVectorFFT should perform a forward FFT. + PS_FFT_REVERSE = 1, ///< psImageFFT/psVectorFFT should perform a reverse FFT. + PS_FFT_REAL_RESULT = 2 ///< 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{verbatim} + +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 and power spectrum must be of the same size as +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{verbatim} +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, psImage *real, const psImage *imag); +psImage *psImageConjugate(psImage *out, const psImage *in); +psImage *psImagePowerSpectrum(psImage *out, const psImage *in); +\end{verbatim} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\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{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{verbatim} +/** A convolution kernel */ +typedef struct { + psImage *image; ///< Kernel data, in the form of an image + int xMin, yMin; ///< Most negative indices + int xMax, yMax; ///< Most positive indices + float **kernel; ///< Pointer to the kernel data + float **p_kernelRows; ///< Pointer to the rows of the kernel data +} psKernel; +\end{verbatim} + +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{verbatim} +psKernel *psKernelAlloc(int xMin, int xMax, int yMin, int yMax); +\end{verbatim} + +\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{verbatim} +psKernel *psKernelGenerate(const psVector *tShifts, const psVector *xShifts, + const psVector *yShifts, bool relative); +\end{verbatim} + +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 a warning shall be generated, following which, +the longer vector trimmed to the length of the shorter, and the +function shall continue. + +\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{verbatim} +psImage *psImageConvolve(psImage *out, const psImage *in, const psKernel *kernel, bool direct); +\end{verbatim} + +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 randon 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{verbatim} +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{verbatim} + +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{verbatim} +psRandom *psRandomAlloc(psRandomType type, unsigned long seed); +void psRandomReset(psRandom *rand, unsigned long seed); +\end{verbatim} + +\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{verbatim} +double psRandomUniform(const psRandom *r); +double psRandomGaussian(const psRandom *r); +double psRandomPoisson(const psRandom *r, double mean); +\end{verbatim} + +\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}. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Astronomy-Specific 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 Metadata +\item Detector and sky positions +\item Astronomy Image +\item Astrometry +\item Photometry +\item Astronomical objects +\end{itemize} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Dates and times} + +\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{timeval} which has been extended to track the time +system being represented. + +\begin{verbatim} +typedef enum { + PS_TIME_TAI, ///< seconds since 1970-01-01T00:00:00Z (Gregorian) + PS_TIME_UTC, ///< seconds since 1970-01-01T00:00:00Z (Gregorian) +} psTimeType; + +typedef struct { + psS64 sec; ///< seconds, negative values represent dates before 1970 + psU32 usec; ///< microseconds + psTimeType type; ///< type of time +} psTime; +\end{verbatim} + +The corresponding constructor shall be: +\begin{verbatim} +psTime *psTimeAlloc(psTimeType type); +\end{verbatim} + +\subsubsection{Current Date and Time} + +Get the current time (in given system): + +\begin{verbatim} +psTime *psTimeGetTime(psTimeType type); +\end{verbatim} + +\subsubsection{Time Conversion} + +Converting between the \code{psTime} time systems is done with: +\begin{verbatim} +psTime *psTimeConvert(psTime *time, psTimeType type); +\end{verbatim} +This function may be used to convert between \code{PS_TIME_TAI} and +\code{PS_TIME_UTC} time representations. The \code{time} is modified +and returned. + +To convert to Local Mean Sidereal Time, it is necessary to provide the +local longitude (specified in radians, positive East of Greenwich) as +well: +% +\begin{verbatim} +psF64 *psTimeToLST(psTime *time, psF64 longitude); +\end{verbatim} +% +The functions may accept either \code{psTimeType}. The \code{time} is +modified and returned. Note that this function must supply the value +UT1-UTC, which is available externally. The value UT1-UTC is +necessary for this an various other SLALIB functions. The following +utility function encapsulates the PSLib mechanism to extract the value +of UT1-UTC: +\begin{verbatim} +double psTimeGetUT1Delta(const psTime *time); +psSphere *psTimeGetPoleCoords(const psTime *time); +\end{verbatim} + +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, +\code{psTimeLeapSeconds}, shall calculate the number of leap seconds +added between two times: +% +\begin{verbatim} +long psTimeLeapSeconds(const psTime *time1, const psTime *time2); +\end{verbatim} + +\subsubsection{External Date and Time Formats} + +A collection of functions convert from the \code{psTime} types to various +external formats. Note that ISO8601 format is "YYYY-MM-DDThh:mm:ss,sZ" + +\begin{verbatim} +psF64 psTimeToJD(const psTime *time); +psF64 psTimeToMJD(const psTime *time); +char *psTimeToISOTime(const psTime *time); +timeval *psTimeToTimeval(const psTime *time); +\end{verbatim} + +\subsubsection{Date and Time Parsing} + +A related collection of functions convert from external formats to a +\code{psTime} type. + +\begin{verbatim} +psTime *psTimeFromJD(psF64 input); +psTime *psTimeFromMJD(psF64 input); +psTime *psTimeFromISO(const char *input); +psTime *psTimeFromTimeval(const timeval *input); +\end{verbatim} + +\subsubsection{Date and Time Math} + +\begin{verbatim} +psTime *psTimeAdd(const psTime *tai1, psF64 delta); +psF64 psTimeSubtract(const psTime *time1, const psTime *time2); +psF64 psTimeDelta(const psTime *time1, const psTime *time2); +\end{verbatim} +\code{psTimeAdd} adds some time to a \code{psTime}. +\code{psTimeSubtract} subtracts two times. \code{psTimeDelta} gives +the absolute time difference between two times. + +Time math is only done on the \code{psTime} TAI type. + +\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 utilise multiple +tables. Some tables may be found at: + +\begin{itemize} +\item \code{ftp://maia.usno.navy.mil/ser7/ser7.dat} +\item \code{ftp://maia.usno.navy.mil/ser7/finals.all} with explanatory + guide at \code{ftp://maia.usno.navy.mil/ser7/readme.finals} . See also + the web page \code{http://maia.usno.navy.mil/}. +\item \code{http://hpiers.obspm.fr/eoppc/eop/eopc01/eopc01.1900-2004} + (contains estimates prior to 1972). +\end{itemize} + +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 MJD, $x_p$ (in arcseconds), $y_p$ +(in arcseconds) and $\Delta$UT (in seconds) separated by whitespace; +blank lines and lines commencing with the comment character (the hash +symbol, \code{#}) as the first character shall be ignored. For +example: + +\begin{verbatim} +# Bulletin A, 9 September 2004. +# ftp://maia.usno.navy.mil/ser7/ser7.dat +# MJD XP(") YP(") UT1-UTC(s) +53258 0.1715 0.4774 -0.45092 +53259 0.1734 0.4757 -0.45039 +53260 0.1753 0.4741 -0.45012 +53261 0.1770 0.4724 -0.45015 +53262 0.1787 0.4708 -0.45045 +\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}). + +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}). + +We define a structure, \code{psTimeTable} to hold a single time table: +\begin{verbatim} +typedef struct { + const char *filename; // Filename of time table + const psF64 validFrom, validTo; // MJDs of validity of time table + bool loaded; // Has the table been loaded? + psVector *mjd; // MJD data, type is psF64 + psVector *dut; // delta UT = UT1 - UTC data + psVector *xp, *yp; // Polar motion data +} psTimeTable; +\end{verbatim} + +The tables shall be read in only when required by the user (hence the +\code{loaded} member); once loaded, the table shall remain in memory +until the termination of the program. The corresponding constructor +shall be: +\begin{verbatim} +psTimeTable *psTimeTableAlloc(const char *filename, psF64 validFrom, psF64 validTo); +\end{verbatim} +which simply constructs a \code{psTimeTable} with the appropriate +\code{filename}, \code{validFrom} and \code{validTo}. Upon +construction, the \code{loaded} member shall be set to \code{false}, +and the vectors containing the data shall be set to \code{NULL}. + +The function \code{psTimeTableLoad} shall load a specified time +\code{table}, returning \code{true} for success, and \code{false} if +the table failed to load. If the \code{table} has already been +loaded, then it shall be re-loaded. +\begin{verbatim} +bool psTimeTableLoad(psTimeTable *table); +\end{verbatim} + +Given a time \code{table} and a Modified Julian Date, \code{mjd}, at +which to interpolate values, \code{psTimeTableInterpolate} shall +calculate the appropriate UT1-UTC (\code{dut}; if non-\code{NULL}), +and polar motion (\code{xp,yp}; if non-\code{NULL}) by interpolation +on the table. If the \code{mjd} is outside the range of the table or +the interpolation was for some other reason unsuccessful, then the +function shall return \code{false}; otherwise the function shall set +those parameters which are non-\code{NULL} (\code{dut,xp,yp}) and +return \code{true}. If the requested \code{mjd} lies in the range of +the \code{table}, but the \code{table} has not been \code{loaded}, +then the function shall call \code{psTimeTableLoad}. +\begin{verbatim} +bool psTimeTableInterpolate(const psTimeTable *table, double *xp, + double *yp, double *dut, psF64 mjd); +\end{verbatim} + +\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} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Regions} + +In many places, we need to refer to a rectangular area. We define a +structure to represent a rectangle: +\begin{verbatim} +typedef struct { + float x0; + float x1; + float y0; + float y1; +} psRegion; +psRegion *psRegionAlloc (float x0, float x1, float y0, float y1); +\end{verbatim} + +\begin{verbatim} +psRegion *psRegionFromString (char *region); +\end{verbatim} +This function converts the IRAF description of a region in the form +\code{[x0,x1:y0,y1]}, used for header entries such as \code{BIASSEC}, +into the corresponding \code{psRegion} structure. + +\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: +\begin{verbatim} +typedef struct { + int id; ///< unique ID for this item + char *name; ///< Name of item + psMetadataType type; ///< type of this item + psElemType ptype; ///< primitive data type + const union { + psS32 S32; ///< integer data + psF32 F32; ///< floating-point data + psF64 F64; ///< double-precision data + void *V; ///< other type + psList *list; ///< psList entry + psMetadata *md; ///< psMetadata entry + } data; ///< value of metadata + char *comment; ///< optional comment ("", not NULL) +} psMetadataItem; +\end{verbatim} + +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{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{psMetadataType}: +% +\begin{verbatim} +typedef enum { ///< type of item.data is: + PS_META_PRIMITIVE, ///< primitive type: use item.ptype + PS_META_LIST, ///< psList; use item.data.list (used for non-unique data) + PS_META_META, ///< psMetadata: use item.data.list + PS_META_STR, ///< string (item.data.V) + PS_META_MATH, ///< psScalar, psVector, psImage (item.data.V) + PS_META_JPEG, ///< JPEG (item.data) + PS_META_PNG, ///< PNG (item.data) + PS_META_ASTROM, ///< astrometric coefficients (item.data) + PS_META_UNKNOWN, ///< other (item.data) + PS_META_NTYPE ///< Number of types; must be last +} psMetadataType; +\end{verbatim} +If the data is a PSLib primitive data value, the primitive data type +is given by the value of \code{ptype}. + +A collection of metadata is represented by the \code{psMetadata} structure: +\begin{verbatim} +typedef struct { + psList *list; ///< list of psMetadataItem + psHash *table; ///< hash table of the same metadata +} psMetadata; +\end{verbatim} +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{psMetadataType} set to +\code{PS_META_ITEM_SET}, and the \code{psMetadataItem.items} element +pointing to a \code{psList} pointing to the actual entries. If +\code{psMetadataItemAlloc} is called with the type set to +\code{PS_META_ITEM_SET}, such a repeated key is created. If the data +value passed to \code{psMetadataItemAlloc} (the quantity in ellipsis) +is \code{NULL}, then an empty \code{psMetadataItem} with the given +keyword is created to hold future entries of that keyword. + +The \code{psMetadataAdd} routine is required to check that all +metadata names are unique unless the type is already qualified as +\code{PS_META_ITEM_SET}; in this case the data are added to the +corresponding \code{psMetadataItem.items} list. + +\subsubsection{Metadata APIs} + +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}-stype formating codes. The +\code{format} entry, which specifies both the metadata type and the +optional flags, is constructed by bit-wise or'ing the appropriate type +and flag. 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{verbatim} +psMetadataItem *psMetadataItemAlloc(const char *name, psMetadataType type, const char *comment, ...); +psMetadataItem *psMetadataItemAllocV(const char *name, psMetadataType type, + const char *comment, va_list list); +\end{verbatim} + +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 using \code{psMetadataItemFree}. +\begin{verbatim} +psMetadata *psMetadataAlloc(void); +\end{verbatim} + +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. Both functions +take an parameter \code{location} which specifies where in the list to +place the element, following the conventions for the \code{psList}. +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_META_ITEM_SET}). +% +\begin{verbatim} +bool psMetadataAddItem(psMetadata *md, psMetadataItem *item, int location); +bool psMetadataAdd(psMetadata *md, int location, const char *name, int format, const char *comment, ...); +\end{verbatim} + +Items may be removed from the metadata by specifying a key or a +location in the list. If the value of \code{name} is \code{NULL}, the +value of \code{location} is used. If the value of \code{name} is not +\code{NULL}, then \code{location} must be set to +\code{PS_LIST_UNKNOWN}. 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. +% +\begin{verbatim} +bool psMetadataRemove(psMetadata *md, int location, const char *key); +\end{verbatim} + +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{verbatim} +psMetadataItem *psMetadataLookup(const psMetadata *md, const char *key); +\end{verbatim} + +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{verbatim} +psMetadataItem *psMetadataGet(const psMetadata *md, int location); +\end{verbatim} + +The metadata may be iterated over by (re-)setting a particular +\code{iterator}, to a location in the \code{psMetadata} list, and +getting the previous or next item. \code{psMetadataGetNext} has the +ability to match the beginning of a key, 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 of metadata --- even those that are non-unique. 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{verbatim} +void psMetadataSetIterator(psMetadata *md, int iterator, int location); +psMetadataItem *psMetadataGetNext(psMetadata *md, const char *match, int iterator); +psMetadataItem *psMetadataGetPrevious(psMetadata *md, const char *match, int iterator); +\end{verbatim} + +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{verbatim} +void psMetadataItemPrint(FILE *fd, const char *format, const psMetadataItem *md); +\end{verbatim} + +\begin{verbatim} +psMetadata *psMetadataReadHeader(psMetadata *output, const char *extname, + int extnum, const char *filename); +\end{verbatim} +Read header data from a FITS image file into a \code{psMetadata} +structure. The \code{extname} and \code{extnum} parameters specify +the extension of interest as above. The data is read from the +extension specified by extname (matching the EXTNAME keyword) or by +the extnum value (with 0 representing the primary header unit (PHU), 1 +the first extension, etc). This function must call \code{psError} and +return \code{NULL} if the specified extension does not exist. If +\code{extname} is not \code{NULL}, \code{extnum} must be set to -1 or +the function shall raise an error. + +\tbd{this function is being deprecated and the equivalent is being +moved to FITS I/O Functions.} + +\begin{verbatim} +psMetadata *psMetadataFReadHeader(psMetadata *output, const char *extname, + int extnum, FILE *f); +\end{verbatim} +Read header data from a FITS image file descriptor into a +\code{psMetadata} structure. + +\tbd{this function is being deprecated and the equivalent is being +moved to FITS I/O Functions.} + +\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{psMetadataParseConfig}, as described below. + +\begin{verbatim} +psMetadata *psMetadataParseConfig(psMetadata *md, int *nFail, const char *filename, bool overwrite); +\end{verbatim} + +Given a metadata container, \code{md}, and the name of a configuration +file, \code{filename}, \code{psMetadataParseConfig} 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}. + +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. + +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}. + +\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}. + +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.} + +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 +by specifying the key at the start of the line (possibly preceded by +whitespace) and specifying the type as a ``multiple symbol'', which we +choose to be an asterisk, \code{*}. No other data may be provided on +this line, though a comment, preceeded 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. + +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: +\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 * # This is a comment in the ITEM_SET item +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. + +We extend \code{psMetadataParseConfig} to allow 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 +interpretted 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{psMetadataParseConfig} 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_META_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}. + +A BNF-like grammar of the configuration file is contained in +\S\ref{sec:configgrammar}. + +\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{verbatim} +typedef struct { + fitsfile fd; +} psFits; +\end{verbatim} +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). + +\subsubsection{FITS Header I/O Functions} + +\begin{verbatim} +psMetadata *psFitsReadHeader(psMetadata *output, const char *extname, int extnum, const char *filename); +psMetadata *psFitsReadHeaderPtr(psMetadata *output, const char *extname, int extnum, const psFits *fd); +\end{verbatim} +Read header data from a FITS image file into a \code{psMetadata} +structure. The \code{extname} and \code{extnum} parameters specify +the extension of interest as above. The data is read from the +extension specified by extname (matching the EXTNAME keyword) or by +the extnum value (with 0 representing the primary header unit (PHU), 1 +the first extension, etc). This function must call \code{psError} and +return \code{NULL} if the specified extension does not exist. If +\code{extname} is not \code{NULL}, \code{extnum} must be set to -1 or +the function shall raise an error. The file is either specified as a +named file on disk or by an open \code{psFits} file pointer. + +\begin{verbatim} +psHash *psFitsReadHeaderSet (const char *filename); +psHash *psFitsReadHeaderSetPtr (const psFits *fd); +\end{verbatim} +Load a complete set of headers from a named file or \code{psFits} file +pointer. This function loads the headers from all extensions into a +\code{psHash}, each entry of which is a pointer to a \code{psMetadata} +structure containing the header data. The hash keys are the +\code{EXTNAME} values for each header (with the value of \code{PHU} +for the primary header unit). + +\begin{verbatim} +bool psFitsWriteHeader(psMetadata *output, const char *filename); +bool psFitsWriteHeaderPtr(psMetadata *output, const psFits *fd); +\end{verbatim} +Write metadata into the header of a FITS image file. The header is +written to the end of the given file. \tbd{Allow for overwriting a +specific header? dangerous in anycase: the header and data must agree +for the file to be valid...} \tbd{consistency check?} \tbd{function +to remove non-FITS header entries?} + +\subsubsection{FITS Image I/O Functions} + +\begin{verbatim} +psImage *psFitsReadImageSection(psImage *output, psRegion region, int z, const char *extname, int extnum, const char *filename); +psImage *psFitsReadImageSectionPtr(psImage *output, psRegion region, int z, const char *extname, int extnum, const psFits *fd); +\end{verbatim} +Read an image or subimage from a named file or \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. +The file pointer version must not make any assumption about the +current position of the file pointer. + +If the native image is a cube, the value of z specifies the requested +slice of the image. The data is read from the extension specified by +extname (matching the EXTNAME keyword) or by the extnum value (with 0 +representing the primary header unit (PHU), 1 the first extension, +etc). 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, if the specified image file does not exist, or 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{verbatim} +bool psFitsWriteImageSection(const psImage *input, int x, int y, int z, const char *extname, int extnum, const char *filename); +bool psFitsWriteImageSectionPtr(const psImage *input, int x, int y, int z, const char *extname, int extnum, const psFits *fd); +\end{verbatim} +Write an image section to the named file (which may exist) or open +\code{psFits} file pointer. This operation may write a portion of an +image over the existing bytes of an existing image. If the file does +not exist, it should be created. If the specified extension does not +exist, it should be created. If an extension is specified and no PHU +exists, a basic PHU should be created. Care must be taken to +interpret x,y,z in the two cases a) there is already an existing image +and b) there is not an existing image. If the image exists, write the +complete psImage data to the existing image starting at the coordinate +x,y,z. If any of these parameters implies writing pixels outside the +existing data area of the image, return an error (ie, if \code{x + +image.nx >= NAXIS1}, \code{y + image.ny >= NAXIS2}, or \code{z >= +NAXIS3}). If the image does not exist, require x,y,z to be zero. +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. + +\tbd{versions of these functions which do not construct the header, +but instead use a corresponding metadata structure, need to exist}. + +\subsubsection{FITS Table I/O Functions} + +\begin{verbatim} +void *psFitsReadTableRow (int *nBytes, char *extname, int row, char *filename); +void *psFitsReadTableRowPtr (int *nBytes, char *extname, int row, psFits *fd); +\end{verbatim} +These functions read a single row of the table in the extension +specified by \code{extname}, in the file specified by the given +\code{filename} or \code{psFits} file pointer. The row number to be +read is given by \code{row}. The result is returned as collection of +\code{nBytes} 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{verbatim} +psArray *psFitsReadTableColumn (char *extname, char *colname, char *filename); +psArray *psFitsReadTableColumnPtr (char *extname, char *colname, psFits *fd); +\end{verbatim} +These functions read a single column of the table in the extension +specified by \code{extname}, in the file specified by the given +\code{filename} or \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 one row's worth of data +per array element. + +\begin{verbatim} +psVector *psFitsReadTableColumnNum (char *extname, char *colname, char *filename); +psVector *psFitsReadTableColumnNumPtr (char *extname, char *colname, psFits *fd); +\end{verbatim} +These functions read a single column of the table in the extension +specified by \code{extname}, in the file specified by the given +\code{filename} or \code{psFits} file pointer. The column is +specified by the FITS table column key given by \code{row}. + +\begin{verbatim} +psArray *psFitsReadTable (int *nBytes, char *extname, char *filename); +psArray *psFitsReadTablePtr (int *nBytes, char *extname, psFits *fd); +\end{verbatim} +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{nBytes}. The result is returned as +a \code{psVector}, which, as a numerical data type, is only valid for +numerical table columns. 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. + +\tbd{need to have corresponding table write functions...} + +\subsection{Detector and Sky 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 for which additional care +must often be taken. We put these into two structures, \code{psPlane} +and \code{psSphere}, respectively. Partitioning these two will enable +error-checking. +% +\begin{verbatim} +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; +\end{verbatim} + +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. + +\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{verbatim} +typedef struct { + psDPolynomial2D *x; + psDPolynomial2D *y; +} psPlaneTransform; +\end{verbatim} + +The \code{psDPolynomial2D} 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] = 0; T.y->mask[1][1] = 0;} + +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{verbatim} +typedef struct { + psDPolynomial4D *x; + psDPolynomial4D *y; +} psPlaneDistort; +\end{verbatim} + +Like \code{psPlaneTransform}, \code{psPlaneDistort} contains two +\code{psDPolynomial4D} 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{verbatim} +psPlaneTransform *psPlaneTransformAlloc(int n1, int n2); +psPlaneDistort *psPlaneDistortAlloc(int n1, int n2, int n3, int n4); +\end{verbatim} + +We require corresponding functions to apply the transformations to a +specified coordinate \code{coords}: +% +\begin{verbatim} +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{verbatim} + +\subsubsection{Celestial Coordinate Conversions} + +We need to be able to convert between ICRS, Galactic and Ecliptic +coordinates, and potentially between arbitrary spherical coordinate +systems. All of these basic spherical transformations represent +rotations of the spherical coordinate reference. We specify a general +transformation function which takes a structure, +\code{psSphereTransform}, defining the transformation between two +spherical coordinate systems (the structure contains the sines and +cosines of the angles involved so as to minimize computation time for +repeated transformations). We also define a function to generate +\code{psSphereTransform}, based on the three angles +describing the location of the pole and the relative equatorial +rotations of the two systems. We also specify special functions to +return the \code{psSphereTransform} for transformations +between standard coordinate systems. + +\begin{verbatim} +typedef struct { + double sinNPlon; ///< sin of North Pole longitude + double cosNPlon; ///< cos of North Pole longitude + double sinNPlat; ///< sin of North Pole lattitude + double cosNPlat; ///< cos of North Pole lattitude + double sinZP; ///< sin of First PT of Ares lon + double cosZP; ///< cos of First PT of Ares lon +} psSphereTransform; +\end{verbatim} + +The constructor is defined as follows: +\begin{verbatim} +psSphereTransform *psSphereTransformAlloc(double NPlon, double NPlat, double ZP); +\end{verbatim} +where \code{NPlon} and \code{NPlat} define the coordinates in the +input system of the north pole in the output system and \code{ZP} +defines the longitude in the input system of the equatorial +intersection between the two systems (e.g, the first point of Ares). +The constructor must calculate the sines and cosines above. + +Spherical coordinates may be transformed by providing the +transformation and the coordinate in the input system to +\code{psSphereTransform}: +\begin{verbatim} +psSphere *psSphereTransformApply(psSphere *out, + const psSphereTransform *transform, + const psSphere *coord); +\end{verbatim} + +The following functions simply return the appropriate +\code{psSphereTransform} to convert between predefined spherical +coordinate systems (i.e., ICRS, Ecliptic and Galactic). These are +constructors as well as the above \code{psSphereTransformAlloc}. +% +\begin{verbatim} +psSphereTransform *psSphereTransformICRSToEcliptic(psTime *time); +psSphereTransform *psSphereTransformEclipticToICRS(psTime *time); +psSphereTransform *psSphereTransformICRSToGalactic(void); +psSphereTransform *psSphereTransformGalacticToICRS(void); +\end{verbatim} + +We also require the ability to precess coordinates from one equinox to +another. + +\begin{verbatim} +psSphere *psSpherePrecess(psSphere *coords, const psTime *fromTime, const psTime *toTime); +\end{verbatim} + +Given coordinates, \code{coords}, with equinox for \code{fromTime}, +the coordinates are precessed to equinox for \code{toTime}. The +\code{coords} are modified in-place. Equinoxes shall be Julian +equinoxes (as opposed to Bessellian). + +\subsubsection{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{verbatim} +typedef struct { + double R, D; ///< coordinates of projection center + double Xs, Ys; ///< plate-scale in X and Y directions + psProjectionType type; ///< projection type +} psProjection; +\end{verbatim} + +The projection type is defined by the following enumerated type \code{psProjectionType}: +\begin{verbatim} +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{verbatim} + +The following functions will project and deproject (respectively) +spherical coordinates: + +\begin{verbatim} +psPlane *psProject(const psSphere *coord, const psProjection *projection); +psSphere *psDeproject(const psPlane *coord, const psProjection *projection); +\end{verbatim} + +\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{verbatim} +psSphere *psSphereGetOffset(const psSphere *position1, + const psSphere *position2, + psSphereOffsetMode mode, + psSphereOffsetUnit unit); + +psSphere *psSphereSetOffset(const psSphere *position, + const psSphere *offset, + psSphereOffsetMode mode, + psSphereOffsetUnit unit); + +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{verbatim} +Note that these should propagate the errors appropriately. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%% Astronomical Images and Astrometry +\include{psLibSDRS_Astrom} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Photometry} + +\tbd{This section is to be deferred, and for now consists only of +placeholders, with no functional items.} + +Photometric observations are performed in an instrumental photometric +system, and must be related to other photometric systems. We +require a data structure which defines a photometric system, as well +as a structure to define the transformation between photometric +systems. + +The photometric system is defined by the psPhotSystem structure. +A photometric system is identified by a human-readable \code{name} +(ie, SDSS.g, Landolt92.B, GPC1.OTA32.r). Each photometric system is +given a unique identifier \code{ID}. Observations taken with a +specific camera, detector, and filter represent their own photometric +system, and it may be necessary to perform transformations between +these systems. Photometric systems associated with observations from +a specific camera/detector/filter combination can be associated with +those components. +\begin{verbatim} +typedef struct { + const int ID; ///< ID number for this photometric system + const char *name; ///< Name of photometric system + const char *camera; ///< Camera for photometric system + const char *filter; ///< Filter used for photometric system + const char *detector; ///< Detector used for photometric system +} psPhotSystem; +\end{verbatim} + +The following structure defines the transformation between two +photometric systems. +\begin{verbatim} +typedef struct { + psPhotSystem src; ///< Source photometric system + psPhotSystem dst; ///< Destination photometric system + psPhotSystem pP, pM; ///< Primary color reference + psPhotSystem sP, sM; ///< Secondary color reference + float pA, sA; ///< Color offset for references + psPolynomial3D transform; ///< Transformation from source to destination +} psPhotTransform; +\end{verbatim} + +The transformation between two photometric systems may depend on the +airmass of the observation and on the colors of the object of +interest. For a specific observation, such a transformations can be +defined as a polynomial function of the color of the star and the +airmass of the observations. If sufficient data exists, the +transformation between the photometric systems may include more than +one color, constraining the curvature of the stellar spectral energy +distributions. This latter term may be significant for stars which +are highly reddened, for example. Derived photometric quantities may +have been corrected for airmass variations, in which case only color +terms may be measurable. The structure defines the transformation +between a source photometric system (\code{src}) and a target +photometric system (\code{dst}). The photometric system of a primary +color is defined by \code{pP, pM} such that the color is constructed +as $pP - pM$. A secondary color is defined by \code{sP, sM}. For +both, a reference color is specified (\code{pA, sA}): the polynomial +transformation terms refer to colors in the form $pP - pM - pA$. The +transformation is specified as a 3D polynomial. For a star of +magnitude $M_{\rm src}$ in the source photometric system, with +additional magnitude information in the other systems $M_{\rm pP}$, +$M_{\rm pM}$, $M_{\rm sP}$, $M_{\rm sM}$, observed at an airmass of +$z$, the magnitude of the star in the target system $M_{\rm dst}$ is +given by: $M_{\rm dst} = M_{\rm src} + transform(z, M_{\rm pP} - +M_{\rm pM} - pA, M_{\rm sP} - M_{\rm sM} - sA)$. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\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} +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 + +\section{Configuration File Grammar} +\label{sec:configgrammar} + +Here follows a BNF-like grammar for the configuration files specified +in \S\ref{sec:configspec}, suitable for use with +\code{Parse::RecDescent}\footnote{\code{http://search.cpan.org/}$\sim$\code{dconway/Parse-RecDescent-1.94/lib/Parse/RecDescent.pod}}. + +\begin{verbatim} +startrule: statement(s) /\z/ + +statement: scalar | vector | multi_declare | comment + +comment: '#' end_of_line + +scalar: name type value comment | name type value + +vector: vname vtype vvalue comment | vname vtype vvalue + +multi_declare: name '*' end_of_line + +name: /[a-z][.\w]*/i + +vname: /\@[a-z][.\w]*/i + +type: + vtype | 'STR' | 'STRING' + +vtype: + 'S8' | + 'S16' | + 'S32' | + 'S64' | + 'U8' | + 'U16' | + 'U32' | + 'U64' | + 'F32' | + 'F64' | + 'C32' | + 'C64' | + 'BOOL' + +value: vector_value | string + +vvalue: vector_value vector_sep vvalue | vector_value vector_sep(?) + +vector_sep: ',' + +vector_value: float | int | bool + +int: integer_constant + +float: floating_constant + +bool: /[tf]\s+?/i + +string: + /\S[^#\n]*/ + +end_of_line: /[^\n]*/ + +### based on syntax found in "C A Reference Manual" + +# integer constants + +integer_constant: + decimal_constant integer_suffix(?) | + octal_constant integer_suffix(?) | + hexadecimal_constant integer_suffix(?) + +decimal_constant: + digit(2..) | + nonzero_digit + +octal_constant: + '0' octal_digit(s) | + '0' + +hexadecimal_constant: + hex_prefix hex_digit_sequence + +digit: + /[0-9]/ + +nonzero_digit: + /[1-9]/ + +octal_digit: + /[0-7]/ + +hex_digit: + /0_9a-f/i + +integer_suffix: + long_suffix unsigned_suffix(?) | + long_long_suffix unsigned_suffix(?) | + unsigned_suffix long_suffix(?) | + unsigned_suffix long_long_suffix(?) + +long_suffix: + /l/i + +long_long_suffix: + /ll/i + +unsigned_suffix: + /u/i + +# floating-point constants + +floating_constant: + decimal_floating_constant | + hexadecimal_floating_constant + +decimal_floating_constant: + digit_sequence exponent floating_suffix(?) | + dotted_digits exponent(?) floating_suffix(?) + +digit_sequence: + digit(s) + +dotted_digits: + digit_sequence '.' digit_sequence | + '.' digit_sequence | + digit_sequence '.' + +exponent: + /e/i sign_part digit_sequence + +sign_part: + /[+-]/ + +floating_suffix: + /[fl]/i + +hexadecimal_floating_constant: + hex_prefix dotted_hex_digits binary_exponent floating_suffix(?) | + hex_prefix hex_digit_sequence binary_exponent floating_suffix(?) + +hex_prefix: + /0x/i + +dotted_hex_digits: + hex_digit_sequence '.' hex_digit_sequence | + '.' hex_digit_sequence | + hex_digit_sequence '.' + +hex_digit_sequence: + hex_digit(s) + +binary_exponent: + /p/i sign_part(?) digit_sequence +\end{verbatim} + +\section{Revision Change Log} +\input{ChangeLogSDRS.tex} + +\end{document}