%%% $Id: pilotADD.tex,v 1.2 2004-04-08 23:49:47 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{00}
\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.

%%%
