Index: /trunk/doc/pilot/ADD/pilotADD.tex
===================================================================
--- /trunk/doc/pilot/ADD/pilotADD.tex	(revision 405)
+++ /trunk/doc/pilot/ADD/pilotADD.tex	(revision 405)
@@ -0,0 +1,862 @@
+%%% $Id: pilotADD.tex,v 1.1 2004-04-08 23:47:32 eugene Exp $
+\documentclass[panstarrs]{panstarrs}
+
+\title{Pan-STARRS IPP Pilot Project}
+\subtitle{Algorithm Definition Document}
+\shorttitle{IPP Pilot ADD}
+\project{Pan-STARRS Image Processing Pipeline}
+\author{Eugene Magnier}
+\group{}
+\organization{}
+\version{DR}
+\docnumber{PSDC-410-002}
+
+\begin{document}
+\maketitle
+
+%% revisions page
+\RevisionsStart
+DR-1 & 2003/12/08 & First draft \\
+\hline
+DR-2 & 2003/12/15 & Second Draft: include various comments, flesh-out algorithms \\
+\hline
+DR-3 & 2003/12/18 & First Release: generally complete, includes reference to spline\_note.pdf \\
+\RevisionsEnd
+
+\DocumentsInternal
+N/A          & Numerical Recipes Press et al \\ \hline
+N/A          & A Note on Splines U.W. AMATH 352 \\ \hline
+NOST 100-1.1 & Definition of the Flexible Image Transport System (FITS) NOST \\
+\DocumentsExternal
+\DocumentsEnd
+
+%% contents page
+\tableofcontents
+\pagebreak
+
+%% figures page
+\listoffigures
+\pagebreak 
+\pagenumbering{arabic}
+
+%%%%% start of main text
+
+\milsection{Scope}
+
+\milsubsection{Identification}
+
+This document establishes the algorithms to be used by the Pan-STARRS
+Pilot Project.  It makes extensive use of algorithms described by
+Numerical Recipes, without sufficient internal references.  
+
+\milsubsection{Document Overview}
+
+Open Issues and TBDs in this document are marked in bold with
+surrounding square brackets.
+
+\milsection{Algorithms}
+
+\milsubsection{Basic Statistics}
+
+Throughput the IPP, we require a number of basic statistical tests and
+numerical methods to be used.  This document defines these algorithims
+in terms of their mathematical operation.  We define the following statistical terms, assuming there is a set of
+data elements $x_i$.
+
+\milsubsubsection{Mean}
+
+The mean is defined as: \[ \bar{x} = \frac{1}{N} \sum_{i = 1}^{N} x_i \]
+
+\milsubsubsection{Clipped Mean}
+
+The clipped mean is used to determine the mean and $\sigma$ of a
+distribution in the presence of outliers. The algorithm for the
+clipped mean has 2 parameters: $N$, the number of iterations and $k$,
+the multiplying factor of the $\sigma$ to exclude outliers.  Typical values
+for $N$ and $k$ are 3 for both.
+
+\begin{enumerate}
+\item Compute the median of the sample. 
+\item Use this median as the first estimator of the mean of the
+  sample: $\bar{x}$
+\item Repeat the following N times:
+\begin{enumerate}
+  \item Compute the standard deviation of the sample $\sigma$.
+  \item Exclude all sample values $x_{i}$ that have $|x_{i} - \bar{x}| > k \sigma$ 
+  \item Compute the mean of the clipped sample
+\end{enumerate}
+\end{enumerate}
+
+\milsubsubsection{Median}
+
+The median is defined as the value for which 50\% of the data values
+are larger and 50\% are smaller.  This can be calculated in two ways
+depending on the size of the data sample.  For small data samples, the
+values are sorted and the median is given by the value of the middle
+element, if the number of values is odd, and by the average of the two
+middle values if the number of values is even.  For large data
+samples, a histogram of values is generated and the bin which contains
+the 50th percentile value is given as the median value.  For integer
+data, the bin size should be an integer.  For floating point data, the
+bin size should be a fraction of the standard deviation: $\sigma / 4$.
+
+\milsubsubsection{Mode}
+\label{mode}
+
+Formally, the mode is the bin in a data distribution with the maximum
+entries.  Care must be taken to define the bin size; given a
+distribution with a single principal component with a distribution
+with standard deviation $\sigma$, a bin size in the vicinity of
+$\sigma$ will properly sample this component of the distribution.  For
+better precision, a bin size a fraction of $\sigma$ may be used and
+the distribution smoothed by the value of $\sigma$.
+
+To measure the mode, determine the value of $\sigma$ for the data
+sample.  Construct a histogram with bins of size $\sigma / 4$ and
+smooth the distribution with a Gaussian of standard deviation $\sigma$
+(see section~\ref{smooth}).  Find the peak bin in the resulting
+histogram: the coordinate of this bin is the mode of the distribution.
+
+\milsubsubsection{Standard Deviation}
+\label{rms}
+
+The standard deviation of the sample is given by:
+
+\[ \sigma = \sqrt{\sum_{i = 1}^N \frac{(x_i - \bar{x})^2}{N - 1}} \]
+
+To minimize the numerical rounding error, this should be calculated
+numerically as:
+
+\[ \sigma = \sqrt{\frac{1}{N - 1} [ \sum_{i = 1}^{N} (x_i - \bar{x})^2 - \frac{1}{N} (\sum_{i = 1}^{N} (x_i - \bar{x}))^2 ]} \]
+
+\milsubsubsection{Chi-Square}
+\label{chisq}
+
+The chi-square is a test of the goodness-of-fit of a model to a set of
+data.  Given a set of $N$ data points $y_i$ with errors $\sigma_i$,
+dependent on free parameters $x_j$, and a model $y (x_j)$ describing
+the relationship between the parameters and the data, the chi-square
+($\chi^2$) is defined as:
+
+\[ \chi^2 = \sum_{i = 1}^N \frac{(y_i - y(x_j))^2}{\sigma^2_i} \]
+
+The reduced chi-square is defined as:
+
+\[ \chi^2_\nu = \frac{\chi^2}{N_{dof}} \]
+where $N_{dof}$, the number of degrees of freedom, is defined as $N -
+N_{par}$ where $N_{par}$ is the number of free parameters $x_j$.
+
+\tbd{Need to point out the purpose of the chi-square as a test of the
+model, and define the confidence limits on the resulting model
+parameters.}
+
+\milsubsection{LU Decomposition}
+\label{LUdecomp}
+
+The LU Decomposition method shall be used to solve matrix equations of
+the form:
+\[ \sum_{j=1}^{N} A_{i,j} \times x_j = B_j \]
+
+We wish to decompose the matrix $A$ with elements $a_{ij}$ into
+diagonal matrices that satisfy the relationship $A = L U$ where $L$ is
+a lower-diagonal matrix of the form:
+\[ L = \left(
+\begin{matrix}
+\alpha_{11} & 0           & 0           & 0 \\
+\alpha_{21} & \alpha_{22} & 0           & 0 \\
+\alpha_{31} & \alpha_{32} & \alpha_{33} & 0 \\
+\alpha_{41} & \alpha_{42} & \alpha_{43} & \alpha_{44} \\
+\end{matrix} \right)
+\]
+%
+and $U$ is an upper-diagonal matrix of the form:
+\[ U = \left(
+\begin{matrix}
+\beta_{11} & \beta_{12} & \beta_{13} & \beta_{14} \\
+0          & \beta_{22} & \beta_{23} & \beta_{24} \\
+0          & 0          & \beta_{33} & \beta_{34} \\
+0          & 0          & 0          & \beta_{44} \\
+\end{matrix} \right)
+\]
+%
+We can find the values of $\alpha_{ij}$ and $\beta_{ij}$ by following
+this procedure.  First, $\alpha_{ii} = 1$.  For all values of $j = 1,
+2, \dots, N$, follow the next steps: For each value of $i = 1, 2,
+\dots, j$, solve for $\beta_{ij}$ using the relationship:
+\[ \beta_{ij} = a_{ij} - \sum_{k=1}^{i-1} \alpha_{ik}\beta_{kj} \]
+For the values of $i = j+1, j+2, \dots, N$, solve for the values of
+$\alpha_{ij}$ with the following:
+%
+\[ \alpha_{ij} = \frac{1}{\beta_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} \alpha_{ik}\beta_{kj} \right) \]
+%
+At this point, use the values of $\alpha_{ij}$ and $\beta_{ij}$ to
+solve for the vector $x_i$ by first solving the equation $L y = B$:
+\[ y_1 = \frac{b_1}{\alpha_{11}}\]
+\[ y_i = \frac{1}{\alpha_{ii}}\left(b_i - \sum_{j=1}^{i-1} \alpha_{ij} y_i\right) \]
+% 
+The values of $y$ may be used to solve for $x_i$ using the
+relationship $U x = b$:
+\[ x_N = \frac{y_N}{\beta_{NN}}\]
+\[ x_i = \frac{1}{\beta_{ii}}\left(y_i - \sum_{j=i+1}^N \beta_{ij} y_ij\right) \]
+
+\milsubsection{Function Minimization}
+
+\milsubsubsection{General Polynomial Fitting}
+
+Given a set of data values $y_i$ with errors $\sigma_i$, related to
+independent data values $x_i$, we would like to fit a polynomial model
+of $N_{par}$ free parameters of the form $y = \sum_{j=0}^{N_{par}}
+\alpha_j x^j$.  The model is determined by minimizing the value of
+$\chi^2$ as defined above (\ref{chisq}), and the minimization is
+determined by solving for the set of parameters $\alpha_j$ for which
+the partial derivatives of the $\chi^2$ with respect to those
+parameters are zero.  The partial derivatives are
+
+\[ \frac{\partial \chi^2}{\partial \alpha_k} = -2 \sum_i (y_i - \sum_j x_i^j \alpha_j) \frac{x_i^k}{\sigma_i^2} \]
+
+Setting the derivatives to zero, this may be reduced to the following
+matrix equation:
+
+\[ \sum_j \alpha_j \sum_i \frac{x_i^k x_i^j}{\sigma_i^2} = \sum_i \frac{x_i^k y_i}{\sigma_i^2} \]
+
+This matrix equation may be solved with LU Decomposition
+(section~\ref{LUdecomp}).
+
+\milsubsubsection{Non-linear Fitting: Levenberg-Marquardt Method}
+
+For models in which the system of equations defined by the partial
+derivatives cannot be solved with the linear technique, other options
+are necessary.  The Levenberg-Marquardt Method (LMM; see Numerical
+Recipes) may be used for these situations.  In LMM, we make a guess at
+the input parameters, measure the $\chi^2$, vary the parameters by a
+particular choice based on the gradient, measure the $\chi^2$ again,
+and adjust the parameters and the parameter varient based on the
+results.
+
+Given a set of $N$ data values $y_i$ with errors $\sigma_i$, dependent
+on values $x_i$, we would like to find the parameters $a_k$ of the
+function $f(x_i; a_k)$ which minimize the $\chi^2$, defined in the
+usual manner (\ref{chisq}).  We start with a set of parameter guesses,
+$a_k$.  We calculate the gradient $\beta_k$ and the Hessian matrix
+$\alpha_{j,k}$ at this parameter selection as follows:
+
+\[ \beta_k = \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{(y_i - f(x_i))}{\sigma_i^2} \] 
+\[ \alpha_{j,k} = \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{\partial f(x_i)}{\partial a_j} \frac{1}{\sigma_i^2} \] 
+%
+We now define the new parameter guess for $a_k$ based on the gradient
+and Hessian by defining $A_{j,k}$ as a variant on $\alpha_{j,k}$ as
+follows:
+
+\[ A_{j,k} = \alpha_{j,k} ~ (j \ne k) \]
+\[ A_{j,k} = (1 + \lambda) \alpha_{j,k} ~ (j = k) \]
+%
+and solve the system of equations represented by:
+\[ A_{j,k} \alpha^\prime_k = \beta_j \]
+%
+where $\alpha^\prime_k$ represents our new attempt at a parameter
+guess. We use this parameter set to calculate $\chi^2$.  If the new
+value of $\chi^2$ is lower than the previous guess, we accept this new
+set of parameters and decrease $\lambda$ by a factor of 10, otherwise
+we keep the old set, and increase the value of $\lambda$ by a factor
+of 10.  We repeat this process until the value of the reduced $\chi^2$
+changes by much less than 1.0.  The resulting values of $a_k$ are the
+best-fit parameters for the system.  If the errors are normally
+distributed, the formal errors on the parameters are then calculated
+by setting $\lambda = 0$ and calculating the covarience matrix
+$C_{i,j}$, the inverse of the matrix $\alpha_{j,k}$.  The independent
+68.3\% confidence limit on parameter $a_k$ is then $\sqrt{C_{k,k}}$.
+Confidence contours for sets of parameters may be defined as well by
+the function $\Delta = \delta\bar{a} P_{j,k}^{-1} \delta\bar{a}$ where
+$P_{j,k}$ is the projected matrix of $C_{j,k}$, ie those rows and
+columns of $C_{j,k}$ associated with the parameters of interest, the
+vector $\delta\bar{a}$.  The value of $\Delta$ is given by the table
+below for specific confidence limits and numbers of parameters.  
+Note that it is necessary to be able to calculate both the function as
+well as its derivative for any combination of parameters and dependent
+variables.
+
+\begin{center}
+\begin{tabular}{|l|r|r|r|}
+\hline
+{\bf P} & \multicolumn{3}{c|}{\bf $N_{par}$} \\
+        & 1    & 2    & 3    \\
+\hline
+68.3\%  & 1.00 & 2.30 & 3.53 \\
+95.4\%  & 4.00 & 6.17 & 8.02 \\
+99.73\% & 9.00 & 11.8 & 14.2 \\
+\hline
+\end{tabular}
+\end{center}
+
+\milsubsection{Smoothing: Boxcar and Gaussian}
+\label{smooth}
+
+Smoothing may occasionally be perfomed on data.  We present the
+algorithms for two typical versions: boxcar and Gaussian smoothing.
+In both smoothing techniques, given a series of data values $f_i$, the
+smoothed values $g_i$ are determined by calculating a linear
+combination based on the input data point and its nearest $2N$
+neighbors in the form:
+
+\[ g_i = \sum_{n=-N}^N c_n f_i \]
+%
+where the values of $c_n$ determine the filter type.  For boxcar
+smoothing, the values $c_n$ are constant, and must be equal to $1/(2N
++ 1)$ to maintain the zeroth moment of the data (care must be taken at
+the ends of the data range to reduce the value of $c_n$ as fewer input
+data points may be used).  For Gaussian smoothing, the crucial
+parameter is $\sigma$, the standard deviation.  The
+value of $N$ should be chosen to be large enough, $N = 5\sigma$, and
+the values of $c_n$ are then just the Gaussian curve:
+
+\[ c_n = \frac{e^{\frac{-n^2}{2\sigma^2}}}{\sqrt{2\pi\sigma^2}} \]
+
+\milsubsection{Robust Statistics Estimations}
+\label{robust}
+
+Given a set of data values $x_i$ with a primary probability
+distribution described by a mean $\bar{x}$ and standard deviation
+$\sigma$ and a much smaller number of outliers with very different
+probability distribution, it is necessary to make estimates of the
+above statistics using a subset of the data to reject the outliers. 
+
+To make a robust measurement of $\sigma$, sort the data values, and
+determine $\sigma_r = (U_{\frac{1}{4}} - L_{\frac{1}{4}}) / 1.34$,
+where $U_{\frac{1}{4}}$ and $L_{\frac{1}{4}}$ are the upper and lower
+quartile values of the distribution (above and below which are found
+25\% of the data points).  Construct a histogram of the data with bins
+of width $\sigma_r / 4$ and smooth the distribution with a Gaussian of
+width $\sigma$.  Find the peak (mode, $\mu$) and fit a Gaussian to the
+histogram portion centered on the mode containing 25\% of the data
+values.  The resulting value $\mu$ is the robust mean and $\sigma$ is
+the robust standard deviation. 
+
+\milsubsection{Splines}
+
+Data interpolation and smoothing may be performed by fitting a spline
+to data points and using the spline function to determine intermediate
+values.  It is necessary to specify the boundary condition choice when
+performing a spline fit.  A very useful and readable formulation for
+three boundary conditions of interest (Natural, Clamped, Not-a-knot)
+is provided by a note from a U. Washington math course note sheet, ``A
+Note On Cubic Splines'' (AMATH 352), included.  \tbd{write this out in
+full later}.
+
+\milsubsection{Data Formats and Conversion}
+
+Image data values are represented in the FITS image in one of several
+possible ways described by the header keyword BITPIX (see FITS
+Standard Document NOST 100-1.1).  The value of BITPIX may be one of:
+
+\begin{center}
+\begin{tabular}{r|r}
+   8 & short integer \\
+  16 & medium integer \\
+  32 & long integer \\
+ -32 & single-precision floating point \\
+ -64 & double-precision floating point \\
+\end{tabular}
+\end{center}
+
+where the absolute value represents the number of bits in the data
+representation and the sign is set to positive for integer
+representations and negative for floating point representations.  The
+byte order is high-significance bytes first (big-endian).  In the
+integer representations, the data array value is represented as a
+signed integer with the number of bits specified by BITPIX.  The
+floating point representations use the IEEE-754 Floating Point
+representation of the array values.  
+
+In addition to the coding described above, the pixel data value may be
+scaled relative to the file data value.  The scaling is determined by
+the FITS header keywords BZERO and BSCALE which define floating point
+coeffiecients of the following equation:
+
+\[ physical = array \times BSCALE + BZERO \] 
+
+where $physical$ represents the actual data value, while $array$
+represents the data value of the bytes in the FITS file.
+
+For the purpose of computation, an internal floating-point
+representation shall be used.  However, data may be converted between
+these representations when written to disk.  Some details of this
+conversion are important.
+
+First, truncation to an integer representation introduces sharp
+transitions in the data which are visible after smoothing.  To avoid
+these artificial transitions, the conversion to an integer value shall
+randomized by 1/2 of that integer value.  Therefore, conversion of a
+physical value to an integer representation using BSCALE and BZERO
+shall use the following formula:
+
+\[ array = (physical - BZERO) / BSCALE + rand(0.5) \]
+
+where $physical$ represents the actual data value, while $array$
+represents the data value of the bytes in the FITS file.
+
+Second, the choice of BZERO and BSCALE shall be set to protect the
+special value of 0.0.  Since the $array$ data values experience
+integer truncation, arbitrary values of BZERO and BSCALE will force an
+initial physical value of 0.0 to have an actual physical value in the
+range -BSCALE / 2 to +BSCALE / 2.  Since 0.0 is a special value,
+BSCALE and BZERO values shall be chosen to force the physical value of
+0.0 to be represented by a number which is in a range of $\pm BSCALE /
+1000$.  This may be achieved by slightly adjusting BSCALE and BZERO
+relative to some initial choice.  The initial choice defines the
+minimum and maximum physical data value ($MIN_p$, $MAX_p$) corresponding
+to the minimum and maximum array data values ($MIN_a$, $MAX_a$) , as
+well as the array data value corresponding to the physical data value of
+0.0 ($ZERO_a$):
+
+\[ MIN_p = MIN_a * BSCALE_i + BZERO_i \]
+\[ MAX_p = MAX_a * BSCALE_i + BZERO_i \]
+\[ 0.0 = ZERO_a * BSCALE_i + BZERO_i \]
+
+For an arbitrary choice of $BSCALE_i$ and $BZERO_i$, the value
+$ZERO_a$ will not be an integer value.  We require the MAX physical data
+value to be maintained, therefore we choose new values of BSCALE and
+BZERO based on the initial values to maintain the value of MAX and to
+force $ZERO_a$ to be the closest integer value:
+
+\[ ZERO_{a,i} = -BZERO_i / BSCALE_i \]
+\[ ZERO_{a,o} = int (-BZERO_i / BSCALE_i) \]
+\[ BSCALE_o = MAX_p / (MAX_a - ZERO_{a,o}) \]
+\[ BZERO_o = - ZERO_{a,o} * MAX_p / (MAX_a - ZERO_{a,o}) \]
+
+With this representation, the physical value corresponding to 0.0 will
+have a value in the range $\pm P$ where $P$ is the precission of the
+keyword representation of BSCALE.  BSCALE shall have a precission of 5
+digits.
+
+\milsubsection{Bias Subtractions}
+
+The subtraction of a bias may be performed in several ways.  This
+document defines the algorithms to be used for these different bias
+subtraction options.  
+
+\milsubsubsection{Overscan Region Definition}
+
+The bias for an image may be determined based on the data values in
+the overscan region.  The overscan region may be defined by the FITS
+header keyword BIASSEC as follows:
+
+\begin{verbatim}
+BIASSEC = '[Xs:Xe,Ys:Ye]'
+\end{verbatim}
+
+where the pixel coordinates are taken with respect to an origin pixel
+at (1,1).  The values {\tt Xs, Xe, Ys, Ye} define a rectangular
+overscan region with {\tt Xs} the first pixel in the x direction, {\tt
+Xe} the last pixel in the x direction, {\tt Ys} the first in the y
+direction and {\tt Ye} the last in the y direction.  Note that the
+single quotes in the definiton of {\tt BIASSEC} are required by the
+FITS standard for string-type keyword values.  
+
+\milsubsubsection{Constant Value Option}
+
+The bias may be represented by a single value determined from the
+entire overscan region pixel values.  These must include:
+
+\begin{itemize}
+\item {\bf mean} of the overscan region
+\item {\bf clipped mean} of the overscan region
+\item {\bf median} of the overscan region
+\end{itemize}
+
+\milsubsubsection{One-Dimensional Function Option}
+
+The overscan region encompases an area which is long in one direction
+compared to the other and which overlaps all data values in that
+direction.  The bias may be represented as a one-dimensional function
+which varies as a function of the pixel coordinate in the long
+direction of the overscan region.  The options for this function
+include:
+
+\begin{itemize}
+\item spline
+\item polynomial
+\end{itemize}
+
+The data values used to determine this function are derived from the
+data in the short dimension of the overscan region, from which a
+single value is derived using one of the estimators:
+
+\begin{itemize}
+\item mean
+\item clipped mean
+\item median
+\end{itemize}
+
+The resulting vector may be further filtered to exclude points which
+had extremely high value of $\sigma$ ($> N \times READ\_NOISE$) or
+which were extreme outliers from the distribution.
+
+\milsubsection{Image Trimming}
+
+A science or calibration image may be trimmed to remove the portions
+of the image which do not contain valid imaging data.  The portion of
+the image which should be retained may be described by the DATASEC
+header keyword as follows:
+
+\begin{verbatim}
+DATASEC = '[Xs:Xe,Ys:Ye]'
+\end{verbatim}
+
+where the pixel coordinates are taken with respect to an origin pixel
+at (1,1).  The values {\tt Xs, Xe, Ys, Ye} define a rectangular
+overscan region with {\tt Xs} the first pixel in the x direction, {\tt
+Xe} the last pixel in the x direction, {\tt Ys} the first in the y
+direction and {\tt Ye} the last in the y direction.  Note that the
+single quotes in the definiton of {\tt DATASEC} are required by the
+FITS standard for string-type keyword values.  
+
+\milsubsection{Flat-field Correction}
+
+The science image should be divided by the flat-field image.  If any
+of the flat-field image pixels have a value less than or equal to zero, the output pixel
+should be set to the bad pixel value, BADPIX.  
+
+\milsubsection{Bad Pixel Masking}
+
+The science image should be left at the original value if the
+equivalent mask pixel value is 0.  Those pixels which have a non-zero
+value in the mask image should be set to the bad pixel value, BADPIX,
+in the science image.
+
+\milsubsection{Sky Subtraction}
+
+The ``sky'', {\it i.e.} the background of the images will be modeled
+for the Pilot project with a two--componenent model: a low spatial
+frequency component, modelled with a 4$^{th}$ order polynomial, and a
+high spatial frequency component corresponding to a fringe image.
+
+Let $(x,y)$ be the coordinates of a pixel in a given cell. The value of
+the background $B(x,y)$ in this pixel is modelized as:
+
+\begin{multline}
+B(x,y)  =   \alpha \times F_{i}(x,y) + C_{0} +C_{1} x + C_{2} y + C_{3} x^2 +
+C_{4} x y + C_{5} y^2 + C_{6} x^3 + C_{7} x^{2} y + C_{8} x y^{y} +
+C_{9} y^3 \\
+        + C_{10} x^4 + C_{11} x^3  y  + C_{12} x^2 y^2 + C_{13} x
+y^{3} +C_{14} y^4 \\ \label{eq:background}
+\end{multline}
+
+where $F_{i}(x,y)$ is the value of the pixel $(x,y)$ in the fringe
+image $F_{i}$. The free parameters of this model are $\alpha$ and
+$C_i$.  There will be a library of $N$ fringe images generated for
+each filter.
+
+In order to recover the background emission and subtract it, the Pilot
+IPP must first make a crude object detection in each cell, as
+described below, in order to perform the sky determination using empty
+sky pixels, then perform a $\chi\text{-squared}$ minimization to
+determine the coefficients of each function.
+
+The algorithm reads as follow:
+
+\begin{enumerate}
+\item make a working copy of the cell image and mask.
+\item measure robust mode and noise ($\sigma$) in the working copy.
+\item mask pixels that have a value greater than the mode plus 5 times
+  the $\sigma$.
+\item using the non masked pixels only, determine the coefficients of
+  the sky model. If there is not enough pixels to perform a sky
+  determination, issue an error.
+\item construct a sky image using the determined coefficients,
+including for those pixels that have been masked in step 3.
+\item subtract this sky image from the original input image.
+\end{enumerate}
+
+Each of the steps are detailed below.
+
+\milsubsubsection{Copy Cell image and Mask}
+
+This is straightforward.
+
+\milsubsubsection{Measure mode and noise.}
+
+The crucial step here is to recognize whether the image is being
+taken on an empty field.
+
+The mode is to be measured using the robust technique describe in
+section~\ref{robust}, using all the non masked pixels of the cell. The
+standard deviation ($\sigma$) is to be measured using the robust
+technique useful in the presence of outliers, described in
+section~\ref{robust}. The $\sigma$ is to be compared to the expected
+$\sigma$, given the mode, the gain of the cell, its readout noise, and
+the subtracted bias value. All these informations are present in the
+metadata. For reminder, the expected $\sigma$ is:
+\begin{equation} 
+\sigma \simeq g \times \sqrt{(O+M/F)/g + (R/g)^{2}}
+\end{equation}
+where $g$ is the detector gain in $e^{-}/DN$, $O$ is the average
+overscan value in DN, M is the mode in DN, F is the flat-field average
+value, $R$ is the readout noise in DN.
+
+Two cases may arise:
+\begin{itemize}
+
+\item if the measured $\sigma$ is more than 3 times the expected
+$\sigma$, this is the sign that the field is crowded (very few pixels
+see the actual sky).  Such an image should be skipped for the pilot
+project.
+\item if the measured $\sigma$ is not more than 3 times the expected
+value, this is the sign that the field is uncrowded and the sky may be
+fitted.
+\end{itemize}
+
+\milsubsubsection{Crude source detection}
+
+Once the mode and $\sigma$ have been determined, a crude object detection
+is performed by checking for each pixel in the cell:
+
+\begin{equation}
+P(x,y) 
+\begin{cases}
+> M + 5 \times \sigma  \longrightarrow \text{Mask pixel}\\
+\le  M + 5 \times \sigma \longrightarrow \text{Do nothing} \\
+\end{cases}
+\end{equation}
+
+where P(x,y) is the value of the image at position (x,y). Note that
+the masking is done on the woking copy of the image. 
+
+\milsubsubsection{Fit the background values}
+
+The fitting method is a standard $\chi-$squared fit. The background B
+is a function of 17 parameters: $B_(\alpha, i, C_{0}, \dots, C_{14})(x,y)$.
+
+The $\chi-$square statistic (see section~\ref{chisq}) with:
+\begin{equation}\label{eq:chisq}
+\chi^{2} 
+ = \sum_{x,y} \left(\frac{P(x,y) - B(x,y)}{\sigma}\right)^{2} 
+ = \frac{1}{\sigma^{2}} \sum_{x,y} \left(P(x,y) - B(x,y)\right)^{2} 
+\end{equation}
+%
+where $\sum_{x,y}$ is a sum over all unmasked pixels in the cell and
+where we use the fact that the $\sigma$ is roughly constant accross the
+cell.
+
+The fit is obtained by finding the  minimum of $\chi^{2}$ in the
+parameter space of $B$. At the minimum, we have:
+\begin{equation}
+\frac{\partial \chi^{2}}{\partial p} = 0
+\end{equation}
+
+where $p$ is one of the 17 parameters. Since the parameter $i$ is a
+discrete parameter (it is the $i^{th}$ fringe frame in the library),
+we cannot compute $\partial \chi^{2} / \partial i$. Instead we loop
+over each fringe image in the library, then form the system of 16
+equations with 16 unknown ($\alpha, C_{0}, \dots, C_{14}$):
+
+\begin{equation}
+\begin{cases}
+\frac{\partial \chi^{2}}{\partial \alpha} = \sum_{x,y} \left(P(x,y) - B(x,y)\right) F_{i}(x,y) = 0 \\
+\frac{\partial \chi^{2}}{\partial C_{0} } = \sum_{x,y} \left(P(x,y) - B(x,y)\right)            = 0 \\
+\frac{\partial \chi^{2}}{\partial C_{1} } = \sum_{x,y} \left(P(x,y) - B(x,y)\right) x          = 0 \\
+\dots \\
+\frac{\partial \chi^{2}}{\partial C_{14}} = \sum_{x,y} \left(P(x,y) - B(x,y)\right) y^4        = 0 \\
+\end{cases}
+\end{equation}
+Explicitly showing the dependence on the parameters, the systems reads:
+
+\begin{equation}
+\setcounter{MaxMatrixCols}{11}
+\begin{matrix}
+\alpha & \sum F_i^2 & + C_0 & \sum F_i & + C_1 & \sum x F_i & + \dots & + C_{14} & \sum y^4 F_i & = & \sum F_i P \\
+\alpha & \sum F_i   & + C_0 &          & + C_1 & \sum x     & + \dots & + C_{14} & \sum y^4     & = & 0 \\
+\alpha & \sum x F_i & + C_0 & \sum x   & + C_1 & \sum x^2   & + \dots & + C_{14} & \sum xy^4    & = & \sum x P \\
+\alpha & \sum y F_i & + C_0 & \sum y   & + C_1 & \sum x y   & + \dots & + C_{14} & \sum  y^5    & = & \sum y P \\
+\vdots &            &\vdots &          &\vdots &            &  \ddots &  \vdots  &              &   & \dots  \\
+\alpha & \sum y^4 F_i & + C_0 &\sum y^4  & + C_1& \sum y^4 x & +\dots & + C_{14} & \sum  y^7    & = & \sum y^4 P \\
+\end{matrix}
+\end{equation}
+
+This system can be solved by LU Decomposition (\ref{LUdecomp}).
+
+Once a solution has been obtained for $(\alpha, C_{0}, C_{1} , \dots ,
+C_{14})$, the $\chi^2$ value is computed, using
+equation~\ref{eq:chisq}. The value of $i$ (i.e. the library template)
+is the one that leads to the lowest $\chi^2$. The goodness of the fit
+should be saved in the metadata. Values of $\chi^2$ that are too large
+are indicative of problems during the fit and warnings must be issued.
+
+If some parameters of the fit are consistant with 0 (e.g. for example
+all degree 3 coefficients compatible with 0), the fit should be
+attempted again without these parameters (e.g.  a polynomial of order
+2 should be fitted), and the goodness of the fit ($\chi^{2}$) should
+be checked against the previous estimate.
+
+\milsubsubsection{Sky image construction}
+
+The critical point here is to reconstruct the sky also behind the
+objects that have been masked in step 3. The sky image $B(x,y)$ is
+reconstructed using equation~\ref{eq:background}, for all the pixels
+in the cell that are not masked in the original input image.
+
+\milsubsubsection{Subtract the sky image}
+
+This is straightforward. Masked pixels in the input image must remain
+masked with the {\tt BADPIX} value.
+
+\milsubsection{Object Detection}
+
+The goal here is to generate a list (catalog) of objects detected as peaks of a
+Gaussian smoothed version of the (assumed to be
+calibrated, masked, and sky-subtracted) source image $f[i_y][i_x]$.
+
+\noindent
+The result shall contain
+
+\begin{itemize}
+\item Locations of the peaks ${\bf r} = (x,y)$ (accurate to sub-pixel
+precision)
+\item The smoothed peak height $f_s$.
+\item An estimate of the `significance' of the peaks $\nu$, defined
+as the ratio of the peak height to the rms value of a smoothed Gaussian white-noise
+image in which the rms pixel value is the same as that estimated for the sky.
+\item The count of masked pixels whose centres lie within a distance $r_{\rm max}$
+of the position ${\bf r}$.
+\item An estimate of the curvature matrix (the second spatial derivatives of
+the smoothed image at the peak).
+\end{itemize}
+
+\noindent
+The input data are:
+
+\begin{itemize}
+\item The rectangular source image $f[i_y][i_x]$ and its dimensions $N_x$, $N_y$.
+\item An associated bad pixel mask $m[i_y][i_x]$.
+\item The dimensions $M_x$, $M_y$ of the zero-padded image on which the FFT operations will be applied.
+\item The Gaussian smoothing scale $r_g$ (in pixels).
+\item The estimated rms sky noise $\sigma$.
+\item The radius $r_{\rm max}$.
+\item The minimum significance threshold $\nu_{\rm min}$.
+\end{itemize}
+
+\noindent
+The main steps are:
+
+\begin{itemize}
+\item Compute the smoothed sky-noise variance $\sigma_s^2$.
+\item Set masked pixels to zero (if not already so).
+\item Pad with zeros.
+\item Compute the complex FFT of the (real) masked and padded image.
+\item Multiply by the Gaussian transfer function $\exp(- k^2 r_g^2 / 2)$.
+\item Inverse transform to obtain the smoothed image $f_s[i_y][i_x]$.
+\item Generate a list of pixels for which $f_s$ exceeds those of the 8 nearest neighbours.
+\item Compute the significance $\nu = f_s / \sigma_s$ .
+\item For each such pixel with $\nu > \nu_{\rm min}$, 
+compute a fractional pixel position correction $d{\bf r}$ using the nearby values.
+\item Compute the count of masked pixels $n_{\rm masked}$ within distance $r_{\rm max}$ of the peak location.
+\end{itemize}
+
+\milsubsubsection{Smoothed Sky Noise Variance}
+
+The smoothed sky-noise rms is to be estimated as
+\begin{equation}
+\sigma_s = \sigma / (2 \sqrt{\pi} r_g)
+\end{equation}
+
+\milsubsubsection{Masking and Padding}
+
+Straightforward.
+
+\milsubsubsection{FFT}
+
+The FFT algorithm should take advantage of the symmetry of the FT of a real
+function: $\tilde f({\bf k}) = \tilde f^\star (-{\bf k})$.
+
+\milsubsubsection{Transfer Function}
+
+The values of the FFT $\tilde f_s$ are defined on a grid with indices $i_x$, $j_x$.
+The spacing of these pixels in $k$-space is $dk_x = 2 \pi / M_x$ and
+$dk_y = 2 \pi / M_y$, so ${\bf k} = (k_x, k_y) = (i_x dk_x, i_y dk_y)$
+from which $k^2 = k_x^2 + k_y^2$.
+
+\milsubsubsection{Fractional Pixel Location}
+
+The integer position ${\bf R}$ is defined to be the location of the
+center of the peak pixel, referred to an origin which is the lower left
+corner of the lower left pixel of the source image.  Note that
+this is not the same as the coordinates in the FITS standard which is perverse,
+and places the spatial origin at (-1/2, -1/2) relative to the origin
+used here.
+
+To obtain the shift in the $x$-coordinate, for a peak pixel with 
+indices $(i, j)$ we define
+\begin{equation}
+F' = (f_s[i+1][j] - f_s[i-1][j]) / 2
+\end{equation}
+and
+\begin{equation}
+F'' = f_s[i+1][j] - 2 f_s[i][j] + f_s[i+1][j]
+\end{equation}
+which are respectively estimates of the gradient and second derivative of the
+smoothed field at the peak.  The corrected position is then
+\begin{equation}
+r_x = R_x - F' / F''
+\end{equation}
+and similarly for the $y$-coordinate.
+Care should be taken to deal with the pathological cases
+where $F''$ is extremely small (in such cases, no sub-pixel
+correction should be applied).
+
+\milsubsubsection{Curvature Matrix}
+
+For a peak at location $(i,j)$ the curvature matrix components are defined as:
+\begin{equation}
+\begin{array}{c}
+f_{sxx} = 2 f_s[i][j] - f_s[i+1][j] - f_s[i-1][j] \cr
+f_{syy} = 2 f_s[i][j] - f_s[i][j+1] - f_s[i][j-1] \cr
+f_{sxy} = f_{syx} = (f_s[i+1][j+1] - f_s[i+1][j-1] - f_s[i-1][j+1] + f_s[i-1][j-1]) / 4
+\end{array}
+\end{equation}
+Together with the peak height $f_s$ and significance these provide estimates of the
+position covariance matrix and also can be used as a descriminator for
+objects which are suitable for astrometry.
+
+% \milsection{Appendices}
+
+\bibliographystyle{plain}
+\bibliography{panstarrs}
+\end{document}
+
+%%%% Spline Description from Robert Lupton %%%%
+Fitting a cubic spline through a set of points ($x_i, y_i$) is well
+known to be a problem that reduces to a tri-diagonal matrix, and
+which can thus be solved in linear time. There are not quite enough
+constraints to make the solution unique, so two extra parameters
+must be specified.
+
+One approach is to specify $d^2y/dx^2$ at the ends of the interval.
+The special choice "$d^2y/dx^2 = 0$" at each end is called a "Natural
+Cubic Spline" and is in general not a good idea, since if $y_i$
+actually have some curvature at the boundaries, then the resulting
+spline approximation will "ring" for the first few gridpoints. E.g. if
+there are 9 points (x,y) in the interval $0,\pi$, then the natural
+spline fit to $\sin(x)$ [with zero-second derivatives at the
+boundaries] has a maximum error of about 6e-5; the fit to $\cos(x)$
+[which has non-zero second derivatives at the boundaries] is about
+8e-3.
+
+A better set of boundary conditions are the so-called "not-a-knot"
+conditions.  In this case, the second derivative is chosen so that the
+spline fitted to the points (i=0, 2, 3, ...) passes through the point
+($x_1, y_1$) [and an equivalent condition at the other boundary].
+
+Another problem with splines is their response to discontinuities in
+the function being fit; the spline function is forced to overshoot,
+and to "ring" for a few points on each side of the boundary. For
+example, a step discontinuity of height 1 at pi/2 with the previous
+9-point sampling produces a maximum overshoot of 0.11.
+
+This can be controlled by the use of `taut splines', which introduces
+extra knots as needed to control such overshooting.
+
+Another useful variation on a cubic spline is a least-squares spline,
+in which case the spline is not required to go through all of the
+points exactly.  In this case, the routine returns a set of knots
+($x_i,y_i$) consistent with the data in a least-squares way.
+
+%%%