%\documentclass[iop,floatfix]{emulateapj}

% \pdfoutput=1

% see latex.readme.txt for notes on using the PS1 template
\documentclass[12pt,preprint]{aastex}
%\documentclass[manuscript]{aastex}
%\documentclass[preprint2]{aastex}
%\documentclass[preprint2,longabstract]{aastex}
\RequirePackage{color}
\input{astro.sty}

% online version may use color, but print version needs b/w
\def\plotmode{col}
%\def\plotmode{bw}

% arxiv needs PDF graphics, but publishers mostly was PS
%\def\plotext{pdf}
\def\plotext{ps}

% use this to make the figure picture path flexible:
%\def\picdir{PATH}
\def\picdir{ALTPATH}

% CZW commands from my previous draft.
\newcommand{\czw}[1]{
  \textbf{CZW: }\textcolor{red}{#1}
}
\newcommand{\czwdraft}[1]{
  \textcolor{red}{#1}
}
\newcommand{\erfcinv}{\mathop{\rm erfcinv}\nolimits}


% Pick a terse version of the title here;
\shorttitle{Pixel Processing in PS1}
\shortauthors{C.Z. Waters et al}
\begin{document}
\title{Pan-STARRS Pixel Processing : Detrending, Warping, Stacking}

% this is a crude trick to get the order of affiliations right.  These
% names are used in the affiliations below.  The user needs to (1) set
% the order and numbers to have the correct sequence in the author
% list and (2) re-order the list at the bottom (and comment-out as needed)
\def\IfA{1}
\def\CfA{2}
\def\MPIA{3}
\def\Princeton{3}
\def\USNO{4}
\def\JHU{1}

% This example has a first author from UH:
\author{
C. Z. Waters,\altaffilmark{\IfA}
IPP Team,
%PS Builder List
% W.~S. Burgett,\altaffilmark{\IfA}
% K.~C. Chambers,\altaffilmark{\IfA} 
% L. Denneau,\altaffilmark{\IfA}
% P. Draper,\altaffilmark{\DUR}
% H.~A. Flewelling,\altaffilmark{\IfA}
% T. Grav,\altaffilmark{\IfA}
% J. N. Heasley,\altaffilmark{\IfA}
% K. W. Hodapp,\altaffilmark{\IfA}
% M. E. Huber,\altaffilmark{\IfA}
% R. Jedicke,\altaffilmark{\IfA}
% N. Kaiser,\altaffilmark{\IfA}
% R.-P. Kudritzki,\altaffilmark{\IfA}
% G. A. Luppino,\altaffilmark{\IfA}
% R. H. Lupton,\altaffilmark{\Princeton}
% E. A. Magnier,\altaffilmark{\IfA}
% N. Metcalfe,\altaffilmark{\DUH}
% D. G. Monet,\altaffilmark{\USNO}
% J.~S. Morgan,\altaffilmark{\IfA}
% P. M. Onaka,\altaffilmark{\IfA}
% P.~A. Price,\altaffilmark{\Princeton}
% C.~W. Stubbs,\altaffilmark{\CfA}
% W.~E. Sweeney,\altaffilmark{\IfA}
% J.~L. Tonry, \altaffilmark{\IfA}
% R. J. Wainscoat,\altaffilmark{\IfA} and 
% C. Z. Waters,\altaffilmark{\IfA}
} % this bracket terminates author list

% The ordering here should be sequential, matching the sequence in the list of authors:
\altaffiltext{\IfA}{Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu HI 96822}
% \altaffiltext{\CfA}{Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138}
% \altaffiltext{\Princeton}{Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA}
% \altaffiltext{\USNO}{US Naval Observatory, Flagstaff Station, Flagstaff, AZ 86001, USA}
% \altaffiltext{\JHU}{Department of Physics and Astronomy, Johns Hopkins University, 3400 North Charles Street, Baltimore, MD 21218, USA}
% \altaffiltext{\MPIA}{Max Planck Institute for Astronomy, K\"onigstuhl 17, D-69117 Heidelberg, Germany}
\begin{abstract}

Lorem ipsum dolor sit amet, consectetur adipiscing elit. Vestibulum
bibendum nisi id tristique posuere. Duis eu mollis nulla. Maecenas est
turpis, mattis tempor urna vitae, placerat rhoncus sem. Lorem ipsum
dolor sit amet, consectetur adipiscing elit. Sed quis velit
nisl. Aliquam erat volutpat. Cras lacinia, nisl tristique auctor
molestie, dolor nulla rhoncus purus, ac accumsan nunc nunc ac
nibh. Maecenas vitae mollis mauris. Ut sollicitudin pulvinar purus,
eget luctus lorem tincidunt vitae. Vestibulum eu mattis neque. Nulla
in tortor id urna dapibus gravida a vel leo.

\end{abstract}

% insert additional keywords as appropriate:
\keywords{Surveys:\PSONE }

\section{OUTLINE}
\begin{verbatim}
* Introduction
* Raw Data Description
* Basic Detrending Steps
  * Overscan
  * Dark
  * Flats
  * Fringes
* Non-traditional / Non-linear issues
  * Persistence & Burntool 
  * faint-end non-linearity
  * regions of bad CTE 
* Variance Maps
* Static Masks
* Dynamic Masks
  * Ghosts
  * Glints
  * Diffraction Spikes
* Magic
* Warping
  * warping kernel
  * linear-by-pieces
  * Covariance 
  * def of skycells?
* Stacking
  * pixel combination rules
  * pixel rejections
  * convolution for matching (success and failure)
* Difference Image analysis
\end{verbatim}

\section{INTRODUCTION}\label{sec:intro}
%% http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?2007ASPC..364..153M&amp;data_type=PDF_HIGH&amp;whole_paper=YES&amp;type=PRINTER&amp;filetype=.pdf
\section{Camera description}

The Pan-STARRS 1 Science Survey uses the 1.4 giga-pixel GPC1 camera with the PS1 telescope on Haleakala Maui to image the sky north of $-30$ declination.  The GPC1 camera is composed of 60 orthogonal transfer array (OTA) devices, each of with is an $8\times{}8$ grid of readout cells.  This parallelizes the readout process, reducing the overhead in each exposure.  However, as a consequence of this large number of individual detector readouts, there are a number of calibrations that need to be included to ensure the response is the same across the entire field of view.

The Pan-STARRS image processing pipeline (IPP) is described elsewhere \citep{MagnierXXX}, but a short summary follows.  The archive of raw exposures is stored on disk, with a database storing the metadata of exposure parameters.  For the PV3 processing, large contiguous regions were defined, and the images for all exposures within that region lauched for the CHIP stage processing.  This stage performs the image detrending (described below in section \ref{dead ref}), as well as the single epoch photometry \citep{MagnierXXY}.  Following the CHIP stage is the CAMERA stage, in which the astrometry and photometry for entire exposure is calibrated against the reference catalog.  This stage also performs masking updates based on the now-known positions and brightnesses of stars that create dynamic features (see \ref{dynamic_masks} below).  The WARP stage is the next to operate on the data, transforming the detector oriented CHIP stage images into sky-oriented images that have common tesselations and sky projections (Section \ref{warping}).  When all WARP stage processing is done for a region of the sky, STACK processing is performed (Section \ref{stacking}) to construct deeper, fully populated images from the set of WARP images that cover that region of the sky.  Beyond the STACK stage, a series of addition stages are done that are described in other papers.  Transient features are identified in the DIFF stage, which takes input WARP and/or STACK data and performs image differencing \citep{HuberXXX}.  Further photometry is performed in the STATICSKY and SKYCAL stages, which add extended source fitting to the point source photometry of objects detected in the STACK images, and calibrate the results against the reference catalog.  The FULLFORCE stage takes the catalog output of the SKYCAL stage, and uses the objects detected in that to perform forced photometry on the individual WARP stage images.  The details of this photometry are provided in \citet{MagnierXXY}.

The full detrend application and processing are described in detail in the the sections below, but a short summary follows.  Once an exposure has been observed on the summit, it is transferred to the main IPP processing cluster at the MRTC-B and registered into the processing database.  This triggers a new chip stage reduction for each of the 60 OTA images that detrends and mosaicks the individual readout cells before measuring the photometric properties of the astronomical objects detected therein.  To begin the detrending, pre-determined static bad pixel masks are used to exclude detector regions that are known to be uncorrectable.  Following this, persisitence trails related to the incomplete transfer of charge in the readout process are corrected.  The image overscan is subtracted, and the known non-linearity of full readouts cells and the substantially worse issues on the edges of the cells are boosted to the expected levels.  The temperature and exposure time dependent dark model is then applied, and the noisemap related to the correlated read noise that is not fully corrected by the dark model is then calculated.  The flat field correction is applied next, and any fringe correction necessary for long wavelength data is subtracted.  Finally, GPC1 specific ``pattern'' corrections are applied to attempt to reduce the cell to cell differences within a single OTA.



% Discuss 2-phase/3-phase device differnces

\section{Burntool / Persistence effect}

Stars that are nearing saturation on GPC1 cause
persistance problems during the read out of the image, creating trails
of light are left on the image.  During the read out process of an
image with a bright star above this threshold, some of the charge
associated with that object is not fully shifted toward the amplifier.
As a result, this charge remains in the starting cell, and is
partially collected in subsequent shifts, resulting in a ``burn
trail'' that extends from the center of the bright source away from
the amplifier (vertically along the pixel columns toward the top of
the cell).

This incomplete charge shifting in nearly full wells continues as each
row is read out.  This results in a remnant charge in the pixels that
the full well was shifted through.  In following exposures, this
remnant charge leaks out, resulting in a trail that extends from the
initial location of the bright source on the previous image towards
the amplifier (vertically down along the pixel column).  This charge
can remain on the detector for up to thirty minutes, requiring the
locations of these ``burns'' needs to be retained between exposures.

Both of these types of persistance trails are corrected via the
BURNTOOL program.  This program does an initial scan of the images,
and identifies stars brighter than a given threshold of 30000 DN.  The
trail from that star is fit with a one-dimensional power law in each pixel column, based on empirical evidence that this
is the functional form of this persistence effect.  Once this fit is
done, the model is subtracted from the image, and the location of the
star is stored in a table along with the exposure PONTIME, which
denotes the number of seconds since the detector was last powered on.

For a subsequent exposure, the table associated with the previous
image is read in, and after correcting trails from the stars on that
new image, it attempts to find remnant trails stored in the table.
These are fit and subtracted using a one-dimensional exponential
model, again based on empirical studies.  If a significant model with
is determined \czwdraft{$\alpha$ < 4}, then this location is retained
in the image output table.  If not, the old burn is allowed to
``expire.''

An issue with this method of correcting the persistance trails is that
it is based on fits to the raw image data, which may have other
signals not determined by the persistence effect.  The presence of
other stars or artifacts along the path of the burn can result in an
incorrect model to be determined, resulting in either an over- or
under-subtraction of the persistance burn. \czwdraft{However, it's
  better than doing nothing.}  

Another issue is that the cores of very bright stars are deformed by
this process, as the burntool fitting preferentially subtracts flux
from one side of the star.  As most stars that result in burns already
have the cores saturated, this does not significantly affect PSF
determination or photometry. \czwdraft{reference to photometry paper?}

\begin{figure}
  \caption{Panel 1: Example image of burn trail.  Panel 2: example image of subsequent image persistence trail.  Panel 3: Repair of panel 1.  Panel 4: Repair of panel 2}
\end{figure}

\begin{figure}
  \caption{example trail data and fit.}
\end{figure}

\section{Masking}
\czwdraft{Technically, we mask the image prior to burntool application now.}

\subsection{Static Masks}

Due to the large size of the detector, it is to be expected that there
will be a number of pixel defects that do not have the detection sensitivity on par
with their neighbors.  To remove these pixels, we have
constructed a static mask that identifies the known defects.  This
mask is constructed in three phases.

First, a CTEMASK is constructed to mask out regions in which the
charge transfer efficiency is low compared to the rest of the
detector.  Twenty-five of the sixty OTAs in GPC1 show some evidence of
CTE issues, with this pattern showing up (to varying degrees) in
roughly triangular patches on the OTA due to defects in the
semiconductor \czwdraft{doping}.  To generate the mask for these
regions, a sample set of \czwdraft{N} evenly illuminated flat field
images were measured to produce a map of the image variance in 20x20
pixel bins.  As the flat image is expected to illuminate the image
uniformly, the expected variances in each bin should be Poissonian
distributed with the flux level.  However, in regions with CTE issues,
adjacent pixels are not independent, allowing the charge in those
pixels to spread.  This reduces the pixel-to-pixel differences,
resulting in a lower-than-expected variance.  All regions with
variance \czwdraft{0.5} smaller than expected are added to the static
CTEMASK.

The next step of mask construction is to examine the flat and dark
models, and exclude pixels that appear to be poorly corrected by these
models.  The darkmask process looks for pixels that are more than
\czwdraft{8} sigma discrepant in \czwdraft{10\%} of the
\czwdraft{test} images after those images have had the dark model
applied to them.  These pixels are assumed to be unstable with respect
to the dark model, and have the DARK bit set in the static mask,
indicating that they are unreliable in scientific observing.
Similarly, the flatmask process looks for pixels that are \czwdraft{3}
sigma discrepant in the same fraction of \czwdraft{test} images after
both the dark and flat models have been applied.  Those pixels that do
not follow the flat field model of the rest of image are assigned the
FLAT mask bit in the static mask, removing the pixels that cannot be
corrected to a linear response.

The final step of mask construction is to examine the detector for
bright columns and other static pixel issues.  This is first done by
processing a set of \czwdraft{100 i filter} science images in the same
fashion as for the darktest.  A median image is constructed from these
inputs along with the per-pixel variance.  These images are used to
identify pixels that have unexpectedly low variation between all
inputs, as well as those that significantly deviate from the global
median value.  Once this initial set of bad pixels is identified, a
$3\times{}3$ pixel triangular kernel is convolved with the initial
set, and any convolved pixel with value greater than \czwdraft{1.0} is
assigned to the static mask.  This does an excellent job of removing
the majority of the problem pixels.  A subsequent manual inspection
allows human interaction to identify other inconsistent pixels
including the vignetted regions around the edge of the detector.
\czwdraft{This might be a lie} As the size of the vignetted region
changes with filter, we have used the g filter to set the baseline
unvignetted region to define the static mask, resulting in the
smallest possible unvignetted region.

\begin{figure}
  \caption{Image map of static mask.  color coded based on mask reason?  It won't be visible at true pixel scale.}
\end{figure}

\subsection{Dynamic masks}

In addition to the static mask that removes the detector level
defects, we also generate a set of dynamic masks that change with the
astronomical features in the image.  These masks are advisory in
nature, and do not completely exclude the pixel from further
processing consideration.  The first of these dynamic masks indicates
the presence of a corrected burntool trail.  These pixels are included
for phtometry, but are rejected more readily in the stacking and
difference image construction, as they are more likely to have small
residual contributions from the under or over subtraction of the
burntool correction.

The remaining dynamic masks are not generated until the IPP camera
stage \czwdraft{IPP paper reference?}, at which point all object
photometry is complete, and an astrometric solution is known for the
exposure.  This added information provides the positions of bright
sources based on the reference catalog, including those that fall
slightly out of the detector field of view or within the inter chip
gaps, where internal photometry may not have identified them.  These
bright sources are the origin for many of the image artifacts that the
dynamic mask identifies and excludes.


\begin{deluxetable}{ccl}
  \tablecolumns{3}
  \tablewidth{0pc}
  \tablecaption{GPC1 Mask Values}
  \tablehead{\colhead{Mask Name} & \colhead{Mask Value} & \colhead{Description}}
  \startdata
  DETECTOR & 0x0001 & A detector defect is present. \\
  FLAT     & 0x0002 & The flat field model does not calibrate the pixel reliably. \\
  DARK     & 0x0004 & The dark model does not calibrate the pixel reliably. \\
  BLANK    & 0x0008 & The pixel does not contain valid data. \\
  CTE      & 0x0010 & The pixel has poor charge transfer efficiency. \\
  SAT      & 0x0020 & The pixel is saturated. \\
  LOW      & 0x0040 & The pixel has a lower value than expected. \\
  SUSPECT  & 0x0080 & The pixel is suspected of being bad. \\
  BURNTOOL & 0x0080 & The pixel may contain an uncorrected or over-corrected burntool streak. \\
  CR       & 0x0100 & A cosmic ray is present. \\
  SPIKE    & 0x0200 & A diffraction spike is present. \\
  GHOST    & 0x0400 & An optical ghost is present. \\
  STREAK   & 0x0800 & A streak is present. \\
  STARCORE & 0x1000 & A bright star core is present. \\
  CONV.BAD & 0x2000 & The pixel is bad after convolution with a bad pixel. \\
  CONV.POOR& 0x4000 & The pixel is poor after convolution with a bad pixel. \\
  MARK     & 0x8000 & An internal flag for temporarily marking a pixel. \\
  \enddata
  \label{tab:mask_values}
\end{deluxetable}
  

\subsubsection{Crosstalk ghosts}

Due to electrical crosstalk between the flex cables connecting the
individual detector devices, ghost objects can be created on some OTAs
due to the presence of a bright source at a different position on the
camera.  Table \ref{tab:crosstalk_rules} summarizes the list of known
crosstalk rules.  In each of these cases, a source object brighter
than -14.47 magnitude (instrumental) creates a ghost object many
orders of magnitude fainter at the target location.  The cell (x,y)
coordinate is identical between source and ghost, as a result of the
transfer occurring as the devices are read.  A circular mask is asdded
to the ghost location with radius $R = 3.44 \left(-14.47 - m_{source,
  instrumental}\right)$ pixels.  Any objects in the photometric
catalog found at the location of the ghost mask have a \czwdraft{flag}
set, marking the object as a likely ghost.  The majority of the
crosstalk rules are bi-directional, with a source in either position
creating a ghost at the corresponding crosstalk target position.  The
two faintest rules are uni-directional, likely due to differences in
the \czwdraft{magical properties of the electronics}.

For the very brightest sources ($m_{instrumental} < -15$), there can
be crosstalk ghosts between all columns of cells during the readout.
These ``bleed'' ghosts were originally identified as ghosts of the
saturation bleeds appearing in the neighboring cells, and as such, the
masking for these objects puts a rectangular mask down from top to
bottom of cells in all columns that are in the same row of cells as
the bright source.  The width of this box is a function of the source
magnitude, with $W = 5 * \left(-15 - m_{source, instrumental}\right)$
  pixels.

\begin{deluxetable}{lllc}
  \tablecolumns{4}
  \tablewidth{0pc}
  \tablecaption{GPC1 Crosstalk Rules}
  \tablehead{\colhead{Type}&\colhead{Source OTA/Cell}&\colhead{Ghost OTA/Cell}&\colhead{$\Delta m$}}
  \startdata
  Inter-OTA & OTA2Y XY3v & OTA3Y XY3v & 6.16 \\
            & OTA3Y XY3v & OTA2Y XY3v &      \\
            & OTA4Y XY3v & OTA5Y XY3v &      \\
            & OTA5Y XY3v & OTA4Y XY3v &      \\
  Intra-OTA & OTA2Y XY5v & OTA2Y XY6v & 7.07 \\
            & OTA2Y XY6v & OTA2Y XY5v &      \\
            & OTA5Y XY5v & OTA5Y XY6v &      \\
            & OTA5Y XY6v & OTA5Y XY5v &      \\
  One-way   & OTA2Y XY7v & OTA3Y XY2v & 7.34 \\
            & OTA5Y XY7v & OTA4Y XY2v &      \\
  \enddata
  \label{tab:crosstalk_rules}
\end{deluxetable}
  
\begin{figure}
  \caption{Figure of crosstalk ghost and bright star source.  Plot of cut across ghost to illustrate the flat-top shape.}
\end{figure}

\subsubsection{Optical ghosts}
% http://arxiv.org/pdf/1207.2513v1.pdf
Due to imperfections in the anti-reflective coating, bright sources
can also result in large out of focus objects, particularly in the
g-filter data.  These objects are the result of light reflecting back
off the surface of the detector, reflecting again off the lower
surfaces of the optics (particularly the L1 corrector lens), and then
back down onto the focal plane.  Due to the extra travel distance, the
resulting source is out of focus and elongated along the radial
direction of the telescope. These optical ghosts can be modeled as a
bright star in location (X,Y) on the focal plane creates a reflection
ghost on the opposite side of the optical axis at (-X,-Y).  The exact
location is fit as a third order polynomial in the focal plane x and y
directions.  An elliptical annulus mask is constructed at the expected
ghost location, with the major and minor axes defined by linear
functions of the ghost distance from the optical axis, and oriented
along the radius of the detector.  All stars brighter than a
filter-dependent threshold (listed in table
\ref{tab:ghost_magnitudes}) have such masks constructed.

\begin{deluxetable}{lc}
  \tablecolumns{2}
  \tablewidth{0pc}
  \tablecaption{Optical Ghost Magnitude Limits}
  \tablehead{\colhead{Filter}&\colhead{$m_{inst}$}}
  \startdata
  g & -16.5 \\
  r & -20.0 \\
  i & -25.0 \\
  z & -25.0 \\
  y & -25.0 \\
  w & -20.0 \\
  \enddata
  \label{tab:ghost_magnitudes}
\end{deluxetable}

\czwdraft{include full polynomial forms?  How best to do that?}

\begin{deluxetable}{lcc}
  \tablecolumns{3}
  \tablewidth{0pc}
  \tablecaption{Optical Ghost Center Transformations}
  \tablehead{\colhead{Polynomial Term}&\colhead{X center}&\colhead{Y center}}
  \startdata 
  $x^0 y^0$ & -1.215661e+02 &  2.422174e+01 \\
  $x^1 y^0$ &  1.321875e-02 &  4.170486e-04 \\
  $x^2 y^0$ & -4.017026e-09 & -1.934260e-08 \\
  $x^3 y^0$ &  1.148288e-10 & -1.173657e-12 \\
  $x^0 y^1$ & -1.908074e-03 &  1.189352e-02 \\
  $x^1 y^1$ &  8.479150e-08 & -9.256748e-08 \\
  $x^2 y^1$ &  1.635732e-11 &  1.140772e-10 \\
  $x^0 y^2$ &  2.625405e-08 &  8.123932e-08 \\
  $x^1 y^2$ &  1.125586e-10 &  1.328378e-11 \\
  $x^0 y^3$ &  2.912432e-12 &  1.170865e-10 \\
  \enddata
  \label{tab:ghost_centers}
\end{deluxetable}

\begin{deluxetable}{lcccc}
  \tablecolumns{5}
  \tablewidth{0pc}
  \tablecaption{Optical Ghost Annulus Axis Length}
  \tablehead{\colhead{Radial Order}&\colhead{Inner Major Axis}&\colhead{Inner Minor Axis}&    \colhead{Outer Major Axis}&\colhead{Outer Minor Axis}}
  \startdata
  $r^0$ & 3.926693e+01 & 5.287548e+01 & 7.928722e+01 & 1.314265e+02 \\
  $r^1$ & 5.325759e-03 &-2.191669e-03 & 1.722181e-02 & -2.627153e-03 \\
  \enddata
  \label{tab:ghost_radii}
\end{deluxetable}

\begin{figure}
  \caption{Figure of full FOV showing optical ghosts.  Possibly only a few OTAs to illustrate shape deformation.}
\end{figure}

\subsubsection{Glints}

\czwdraft{I thought we stopped this because of a hardware change?  Is
  that not true?}  Prior to \czwdraft{DATE}, a reflective surface at
the edge of the camera aperture was open to light passing through the
telescope.  Sources brighter than $m = -20$ that fell on this
reflective surface resulted in light being scattered across the
detector surface in a long narrow glint.  This surface was physically
masked on \czwdraft{DATE} \czwdraft{right?}, but data prior to that
have a dynamic mask constructed when a reference source falls on the
focal plane within \czwdraft{approximately} one degree of the detector
edge.  This mask is 150 pixels wide, with length $L = 2500 \left(-20 -
m_{inst}\right)$.  \czwdraft{Am I correct that this is basically a one-degree edge around the detector?}

%%
%% GLINT_MAX_MAG                   F32 -21.0
%% GLINT.REGION                    MULTI

%% GLINT.REGION                    METADATA
%%   REGION                        STR  [-38000:-24000,-20000:+20000]
%%   GLINT.TYPE                    STR  LEFT
%% END

%% GLINT.REGION                    METADATA
%%   REGION                        STR  [+24000:+38000,-20000:+20000]
%%   GLINT.TYPE                    STR  RIGHT
%% END

%% GLINT.REGION                    METADATA
%%   REGION                        STR  [-20000:+20000,+24000:+38000:]
%%   GLINT.TYPE                    STR  TOP
%% END

%% GLINT.REGION                    METADATA
%%   REGION                        STR  [-20000:+20000,-38000:-24000]
%%   GLINT.TYPE                    STR  BOTTOM
%% END

\begin{figure}
  \caption{Example of glint.}
\end{figure}

\subsubsection{Diffraction spikes}

Bright objects also form diffraction spikes that are dynamically
masked.  These are filter independent, and are modelled as rectangles
with length $L = 10^{0.096 * (7.35 - m)} - 200$ and width $W = 8 + (L
- 200) * 0.01$.  These spikes are dependent on the camera rotation,
and are oriented at $\theta = n * \frac{\pi}{2} - \mathrm{ROTANGLE} +
0.798$.

\subsubsection{Saturated stars}

The cores of saturated stars are masked as well, with radius $r = 10.15 * (-15 - m_{inst})$.  \czwdraft{good job here.}

\begin{figure}
  \caption{Example of saturated star, which will also nicely show the diffraction spikes.}
\end{figure}

\subsection{Video Mask}

One aspect of the OTAs in GPC1 is that an individual cell can be read
off repeatedly while the other cells integrate, resulting in a video
signal from that cell.  This data is used for telescope guiding
purposes, and a single exposure is likely to have a number of these
video cells in different OTAs.  However, reading these cells while
integrating on the others changes the characteristic dark model (see
below) experienced by the other cells on the OTA.  The observed effect
of this is that the glow associated with the amplifiers in the corners
of the cells is depressed during the video readout, relative to the
nominal glow.  Because of this, the standard dark model oversubtracts
this glow.  Before the nature of this issue was fully understood,
these poorly constrained corners were masked with 25-pixel radius
quarter circles, centered on the (0,0) pixel nearest the cell
amplifier.  The other corners of the cell were masked with a 15-pixel
radius quarter circle, as the amplifier location is off the edge of
the cell.


\subsection{Masking fraction}

For the full field of view that falls on the sixty OTAs, 14.7\% \czwdraft{check this} of all pixels are masked.  The majority of this masking is due to regions that fall within the vignetted region.  Defining the radius of the unvignetted region to be 3 degrees, and excluding pixels that fall beyond this point reduces the static masking fraction to 9.7\%.

Unfortunately, due to the design of the OTAs and readout cells, a non-negligible fraction of the field of view falls onto an area that does not have a detector pixel.  For a given OTA mosaicked to a $4846\times{}4868$ pixel image, the 64 $590\times{}598$ pixel readout cells cover 95.7\% of the OTA area.  

For the inter-chip gap area loss, we use two field of view calculations.  The reference field of view of GPC1 is 3 degrees, which at the nominal plate scale of 0.258 arcseconds per pixel, translates to a 20930 FPA pixel radius.  However, based on the manual masking of the vignetted region, illuminated pixels are generally unvignetted out to 3.25 degrees, or a 22720 FPA pixel radius.  Although these result in different coverage areas, summing the number of pixels on OTA mosaicked images ($4846\times{}4868$ pixels) within either field of view results in a inter-chip gap mask fraction of 7\%.

%% mysql> select filter,AVG(camProcessedExp.maskfrac_ref_static), AVG(camProcessedExp.maskfrac_ref_dynamic), AVG(camProcessedExp.maskfrac_ref_advisory), AVG(camProcessedExp.maskfrac_max_static),AVG(camProcessedExp.maskfrac_max_dynamic),AVG(camProcessedExp.maskfrac_max_advisory) from camRun join camProcessedExp USING(cam_id) JOIN chipRun USING(chip_id) JOIN rawExp USING(exp_id) WHERE camRun.label = 'LAP.PV3.20140730.final' GROUP BY filter;
%% +---------+------------------------------------------+-------------------------------------------+--------------------------------------------+------------------------------------------+-------------------------------------------+--------------------------------------------+
%% | filter  | AVG(camProcessedExp.maskfrac_ref_static) | AVG(camProcessedExp.maskfrac_ref_dynamic) | AVG(camProcessedExp.maskfrac_ref_advisory) | AVG(camProcessedExp.maskfrac_max_static) | AVG(camProcessedExp.maskfrac_max_dynamic) | AVG(camProcessedExp.maskfrac_max_advisory) |
%% +---------+------------------------------------------+-------------------------------------------+--------------------------------------------+------------------------------------------+-------------------------------------------+--------------------------------------------+
%%             static              dynamic                advisory
%% | g.00000 |   0.19642137972007 | 0.00010322263512709 |    0.026838445469766 
%%           |   0.20949461794863 |   9.89200027293e-05 |    0.026431927734548 | 
%% | r.00000 |   0.19675996201399 | 0.00025214447869606 |    0.032641054600788 
%%           |   0.20989768279138 | 0.00023994155711801 |    0.032178525485201 | 
%% | i.00000 |   0.19677587604327 | 0.00057470697316504 |    0.038096251937072 
%%           |   0.21003570722292 | 0.00053987093278142 |    0.037471018638997 | 
%% | z.00000 |    0.1974290315691 | 0.00024758901226967 |     0.03064123748973 
%%           |   0.21055007930696 | 0.00023452690039757 |    0.030144453360769 | 
%% | y.00000 |   0.19828990634315 | 0.00014523787521897 |    0.021984846417987 
%%           |   0.21130344126869 | 0.00013634812877977 |     0.02163070300815 | 

Summing mask fractions from these three contributions results in an average of $\sim 20\%$ masking fraction across the field of view.  Dynamic masking adds an additional $2-3\%$, with advisory burntool masking contributing the largest component.

\section{Overscan}

Each cell on GPC1 has an overscan region that covers the
first\czwdraft{?} 34 columns of each row, and the last\czwdraft{?} 10 rows
of each column.  No light lands on these pixels, so the image region
is trimmed to exclude them.  Each row has an overscan value
subtracted, calculated by finding the median value of that row's
overscan pixels.  These medians are then smoothed between rows with a
3-row wide boxcar.  

\section{Non-linearity Correction}

The pixels of GPC1 are not perfectly linear at all flux levels.
Particularly, at low flux levels, some pixels have a tendency to sag
relative to the expected linear value.  This effect is most pronounced
along the edges of the detector, although some entire cells show
evidence of this effect.

To correct this sag, we study the flux behavior of a series of dark
frames with a ramp of exposure times.  As the exposure time increases,
the flux on each pixel also increases in what is expected to be a
linear manner.  Each of these dark exposures in this exposure time ramp is overscan
corrected, and then the median is calculated for each cell, as well as for 
the rows and columns within ten pixels of the edge of the science
region.  From these median values at each exposure time value, we can
construct the expected trend by fitting a linear model, $f_{region} =
gain * t_{exp} + bias_0$, to the median fluxes for darks with exposure
times between 3 and 12 seconds.  This time interval was selected as it
avoids the non-linearity at low fluxes, as well as the possibility of
high-flux non-linearity effects.  From this set of models for each
row, column, or full cell, we construct a table of correction values
by linear interpolating the row and column results to match the full
cell results in the center of the detector.

This non-linearity effect appears to be stable in time, with little
evident change over the survey duration.

\czwdraft{I have figures at http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearity that might be useful}
%http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearity_AllEdges
%http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/DetectorLinearityArchive

\begin{figure}
  \caption{Example plot of linearity as a function of incident brightness.}
\end{figure}

\section{Dark/Bias Subtraction}
% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/Background_Dark_Model
The dark model we make for GPC1 considers each pixel individually,
independent of any neighbors.  To create the dark model, we fit an multi-dimensional model to the array of input pixels
from a randomly selected set of 100-150 \czwdraft{overscan corrected}
dark frames chosen from a given date range.  The model fits
each pixel as a function of the exposure time $t_{exposure}$ and the
detector temperature $T_{chip}$ such that $dark = a_0 + a_1
t_{exposure} + a_2 T_{chip} t_{exposure} + a_3 T_{chip}^2
t_{exposure}$.  This fitting is performed over the sample of input pixels,
and the coefficients $a_i$ stored in the detrend image.  The constant
$a_0$ term includes the bias signal, and as such, a separate bias
subtraction is not necessary.

Applying the dark model is simply a matter of calculating the response
for the exposure time and detector temperature for the image to be
corrected, and subtracting the resulting dark signal from the image.

\subsection{Time evolution}

\czwdraft{The dark model is noticably unstable on time scales of months, and so we have generated a sequence in time to keep the effect of a missed correction low.}

The dark model is not consistently stable over the full survey, with significant drift over the course of multiple months.  Some of the changes in the dark can be
attributed to changes in the voltage settings of GPC1, but the
majority seem to be the result of some unknown parameter.  We
can separate the dark model history of GPC1 into three epochs.  The
first epoch covers all data taken prior to 2010-01-23.  This epoch
used a different header keyword for the detector temperature, making
data from this epoch incompatible with later dark models.  

The second epoch covers data between 2010-01-23 and 2011-05-01, and is
characterized by a largely stable but oscillatory dark solution.
There appear to be two modes that the dark model switches between
apparently at random.  No clear cause has been established for the
switching, but there are clear differences between the two modes
\czwdraft{figures?}.

The evidence of these two modes comes from the discovery of a slight
gradient along the rows of certain cells.  This is a result of a drift
in the bias level of the detector.  Therefore, an appropriate dark
model should remove this gradient entirely.  For these two modes, the
magnitude of this bias drift is different, so a single dark model over
corrects the low-magnitude mode, and undercorrects the high-magnitude
mode.  Upon identifying this two-mode behavior, and determining the
switching points, two separate darks models were constructed from
appropriate ``A'' and ``B'' mode dark frames.  Using the appropriate
dark minimizes the effect of this bias gradient in the dark corrected
data.  Table \ref{tab:dark mode dates} lists the dates used for each dark mode.

\czwdraft{The evidence of the mode switching can be visualized in Figure \ref{fig:dark switching}.  This figure shows image profile along the x-pixel axis binned along the full y-axis of dark corrected images for OTA67.  These images are from sequential days, and have been corrected with a dark model constructed from the full set of dark data within the second epoch.  The opposite sign of the slopes of these profiles indicates that the average dark model does not correct these dates sufficiently, due to the contradictory dark signals between the two modes.}

After 2011-05-01, the two-mode behavior of the dark disappears, and is
replaced with a slow dateobs-dependent drift in the magnitude of the
gradient.  This drift is sufficiently slow that we have modeled it
using three dateobs-independent dark model for different date ranges.
These darks cover the range from 2011-05-01 to 2011-08-01, 2011-08-01
to 2011-11-01, and 2011-11-01 and on.  The reason for this time
evolution is unknown, but we seem to be able to model it with
reasonable accuracy by creating new dark models.

\begin{figure}
  \caption{Example of raw and dark calibrated exposure.  Plots of horizontal cuts for A/B/average corrections.}
\end{figure}

\subsection{Video Dark}

Dark signal is stronger in cell corners due to amplifier glow.  Standard model corrects this.  When OTA reads video cell, the dark model changes.  The standard model oversubtracts the dark model.  Make video darks from dark data that has had video signal running.  Need two passes to cover all cells (shifting video cell between the two).  Can construct the video dark and the standard dark simultaneously, by using OTAs that have video on and off.

Video dark signal appears linear and stable, allowing archival data from prior to video dark data to be corrected by simply taking $VD_{2009} = D_{2009} - D_{Modern} + VD{Modern}$.

\section{Noisemap}

Based on a study of the positional dependence of detected sources, we have discovered that the cells in GPC1 do not have uniform noise characteristics.  Instead, there is a gradient along the pixel rows, with the noise generally higher away from the read out amplifier.  This is likely another effect of the row-by-row bias issue discussed below.  This gradient has the effect that the read noise increases as the row is read out.  As a result of this increased noise, more sources are detected when the readnoise is assumed constant across the readout.  To mitigate this noise gradient, we construct a set of noisemap images by measuring the median variance on bias frames.  The variance is calculated in boxes of 20x20 pixels, and then linearly interpolated to cover the full image.  

Unfortunately, due to correlations in the row-to-row offsets \czwdraft{in the noise?}, the variance measured from the bias images does not fully remove the positional dependence of objects that are detected.  The reason for this is that the simple noisemap underestimates the noise observed when the image is filtered during the object detection process.  This filtering convolves the background noise with a PSF, which has the effect of amplifying the correlated peaks in the noise.  This amplification can therefore boost background fluctuations above the threshold used to select real objects, contaminating the final object catalogs.

To resolve this issue, we chose a PSF with a typical FWHM, and used it to look for detections on a sample of bias images.  As the bias has no real sources, all objects found are by definition false, which provides an idea of how much our noisemap estimation deviates from the ``true'' noise observed by the object detection process.  For a region of area $X*Y$, if we find $k$ false detections above our signal-to-noise threshold $sigma_{thresh}$, then we can estimate how much the noise model deviates from what is observed.  The observed noise threshold is defined as $\sigma_{observed} = \sqrt{2} * \erfcinv{2 * k A_{psf} / (X * Y * N_{exp})}$, where $A_{psf}$ is the footprint size of the PSF (taken as 16 pixels), and $N_{exp}$ is the number of exposures examined in this location.  From this observed threshold, we scale the noisemap previously calculated by the boost factor $B = \sigma_{thresh} / \sigma_{observed}$.  

The row-to-row variations that contribute to the extra noise are related to the dark model, and because of this, as the dark model changes, the effective noise also changes.  To ensure that the noisemap accurately matches the true noise level, we have created different noisemap models for the three major time ranges of the dark model.  We do not see any evidence that the noisemaps have the A/B modes visible in the dark, and so we do not generate different models for each individual dark model.  

\section{Remnance?}

\czwdraft{Despite the known persistence effects of the detectors, we do not do any remnance correction beyond what is discussed above in the burntool section.  Therefore, I probably should just remove this section entirely.}

\section{Shutter?}

\czwdraft{I don't believe that we do a shutter correction either.  So, again, probably shouldn't include it.}

\section{Flat}

\czwdraft{I don't know how the flat calibration code works.  We start with flat field images of the sky, but due to the size of the detector, it is difficult to equally illuminate each pixel.  Therefore, flat calibration.}

Determining a flat field correction for GPC1 is a challenging
endeavor, as the wide field of view makes it difficult to construct a
uniformly illuminated image.  Using a dome screen is not possible, as
the variations in illumination and screen rigidity create unusably
large scatter between different images that are caused by the detector response function.  Because of this, we use sky
flat images taken at twilight, which are more consistently illuminated
than screen flats.  We calculate the mean of these images to determine
the starting flat model.

From this initial flat model, we construct a correction to remove the
effect of the problems illuminating the large area.  This is done by
dithering a series of science exposures across a given pointing.  By
comparing the measured fluxes for a given star as a function of
position on the detector, we can determine the position dependent
scaling factors.  These scale factors can then be used to correct the
initial flat field model to better represent the detector response.

The flat model appears stable with time, although directly measuring
this is as difficult as originally constructing the model.  However,
due to the photometric consistency observed in the catalog of GPC1 measurements, we
can be confident that the flat model is not as time dependent as the
dark correction.


\section{Pattern correction}

Due to the row-by-row bias offsets that are not cleanly removed by the
dark model, we have a set of ``pattern'' corrections that are applied
to some selection of the images.  The PATTERN.ROW correction is used
to remove the remaining row-by-row variation, and the PATTERN.CELL and
PATTERN.CONTINUITY corrections attempt to ensure that the cells of a
given OTA are consistent with the other cells on that OTA.  These corrections are
largely designed to fix issues that are not stable enough with time
for the dark model or flat field model to fully account for the
detector behavior.

\subsection{Pattern Row}
% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/GPC1_Bias_Pattern_Study
As discussed above in the dark and noisemap sections, certain
detectors have significant row-by-row bias offsets.  The magnitude of
these offsets increases as the distance from the readout amplifier
increases, resulting in horizontal streaks that are more pronounced
along one edge of the cell.  As the level of the offset is largely
random, the dark correction cannot fully remove this structure from
the images, and the noisemap level only indicates the level of the
variance added by these bias offsets.  Therefore, we apply the
PATTERN.ROW correction in an attempt to mitigate the offsets.  To
force the rows to agree, a \czwdraft{second} order polynomial is fit to
each row in the cell, and that trend subtracted from the data.  The
median offset (corresponding to the background level) for each row is then fit by a first order polynomial, and that trend is then added
back to the image so that the sky level on the cell matches its neighbors during
background subtraction.

This correction was required on all cells on all OTAs prior to
\czwdraft{2009-12-01}, at which point a modification of the camera
electronics resolved the row-by-row offsets for the majority of the
detectors.  As a result, we only apply this correction where it is
necessary, as shown in Figure \ref{fig: pattern row required}.

Although this correction does resolve the row-by-row offset issue in a
satifactory way, large and bright astronomical objects can bias the
fit significantly.  This results in an oversubtraction of the offset
near these objects.  As the offsets are calculated on the pixel rows,
this oversubtraction is not uniform around the object, but is
preferentially along the horizontal x axis of the object.  

%% \czwdraft{keep this?}  This row-by-row offset is visible in similar
%% camera designs, and has been removed by identifying the noise signal
%% in the pixel data stream.  By taking the FFT of the pixels and a
%% reference signal, the frequency of this noise can be isolated and
%% removed, resulting in a much cleaner image.  However, GPC1 does not
%% record the value of the reference signal, instead automatically
%% subtracting it from the data values.  Without this comparison signal,
%% we have been unable to reproduce this method, as there is no obvious
%% FFT component visible.

\begin{figure}
  \caption{Example of pre/post pattern row application.}
\end{figure}

\subsection{Pattern Cell}

As the bias level of a given cell may not exactly match that of its
neighbors, fitting a smooth background model results in over and
under-subtraction of the sky level at the cell boundary
discontinuities.  The PATTERN.CELL correction was the first attempt to
remove this effect on the worst cells, by forcing all the cells of an
OTA to the same level.  Each cell has the median value measured, and
then each cell has an offset added that shifts the cell to match the
median of those medians.

This correction is reasonable when the astronomical signal is smooth,
with no objects that are large relative to the size of an individual
cell.  However, the presence of large galaxies (or even bright stars)
can bias the offsets for some cells from their neighbors.  Because of
this issue, we no longer apply this correction to any data.

\subsection{Pattern Continuity}

As the PATTERN.CELL correction was clearly insufficient in many
situations, we designed a replacement correction that would lower the
distortion for large objects less.  In addition, studies of the
background level illustrated that the row-by-row bias introduces
small background gradient variations along the rows of the cells that is not stable enough to be completely fit by the dark model.  This results
in a ``sawtooth'' pattern horizontally across an OTA, and as the background model
assumes a smooth sky level, this induces over and under
subtraction at cell boundaries.  As the PATTERN.CELL was designed to
correct mean changes between cells, it could not adequately resolve
this higher order issue.

The replacment for PATTERN.CELL was the PATTERN.CONTINUITY correction,
which attempts to match the edges of a cell to those of its neighbors.
For each cell, a thin box \czwdraft{10} pixels wide on each edge is extracted and the median
value calculated for that box.  These median values are then used to
construct a vector of differences $diff_i = \sum_{j,j'} Edge_{i,j} -
Edge)_{i',j'}$, along with a matrix of associations $A_{i,i'} =
\sum_{j,j'} \delta(j,j')$ denoting which cell boundaries touch
another.  By solving the system $A x = diff$, we can find the set of
offsets $x_i$ that should be applied to each cell to ensure the
minimum differences between cells.

Due to the known slope in some cells, the effect of this correction is
to align the cells into a single ramp, at the expense of the absolute
background level.  However, as we subtract off a smooth background
model, the deviations from an absolute sky level are unimportant.  The fact that the final
ramp is smoother than it would be otherwise also allows for the
background subtracted image to more closely match the astronomical
sky, without over- and under-subtractions at cell edges.

\begin{figure}
  \caption{Continuity example, with background issue.}
\end{figure}

\section{Fringe correction}

Due to variations in the thickness of the detectors, we observe
interference patterns at the infrared end of the filter set, as
the wavelength of the light becomes comparable to the thickness of these variations.
Visually inspecting the images shows that the fringing is most
prevalent in the y-filter images, with minimal fringing in other
bands.  As a result of this, we only apply a fringe correction to the y filter
data.

The fringe is constructed by randomly determining a set of boxes for
each OTA cell, and measuring the sky subtracted median value in those
boxes for a series of images.  These samples are selected at the same
location on each image, allowing the astronomical signal to be
filtered out as an additional noise term.  A least squares fit to the
data is then calculated, providing the model of the fringe strength at
that location.

Applying the fringe is done in the same way, with samples measured
across the image to determine the relative strength of the fringing in
this image.  The solution derived from the detrend is then scaled to
match that observed in the science image, and subtracted away.

\section{Background subtraction}

\czwdraft{A background model is generated for each OTA, once all the individual cells have been mosaicked together.  Super-pixels are then defined that divide the image into XxY subregions, and the mean calculated for each subregion.  This grid is shifted by a half-width, and the means recalculated, to double the sampling frequency.  A background model is then calculated by interpolating over this sampled grid.}

%% * Magic
%% * Warping
%%   * warping kernel
%%   * linear-by-pieces
%%   * Covariance 
%%   * def of skycells?
%% * Stacking
%%   * pixel combination rules
%%   * pixel rejections
%%   * convolution for matching (success and failure)
%% * Difference Image analysis


\section{Warping}

To provide a consistent and uniform set of images for co-added image
stacking and image differences, the individual mosaicked OTA images
are projected onto a common set of tangent plane projected regions.
These projection cells are $4\times{}4$ degree fields spaced onto set
of projection centers that fully cover the sky.  These projection
cells are arranged into rings of constant declination, and allowed to
overlap as $|\delta|$ increases.  Each projection cell is further
subdivided into \czwdraft{size} sky cells, which have constant overlap
regions of \czwdraft{overlap}.  These skycells are the main image unit
used for processing image data beyond the initial chip stage.

After the detrending and photometry, the detection catalog for the
full camera is fit to the reference catalog, producing third-order
astrometric solutions that map the detector focal plane to the sky,
and map the individual OTA pixels to the detector focal plane.  This
solution is then used to determine which skycells the exposure OTAs
overlap.

Foreach output skycell, all overlapping OTAs and the calibrated
catalog are read into the \textbf{pswarp} program.  Each input image
is examined in order, and the same transformation performed.  This
transformation breaks the output warp image into $128\times{}128$
pixel grid boxes.  Each grid box has a locally linear map calculated
that converts the output warp image coordinates to the input chip
image coordinates.  By doing the transformation in this direction,
each output pixel has a unique sampling position on the input image
(although it may be off the image frame and therefore not populated),
preventing gaps in the output image due to the spacing of the input
pixels.

With the locally linear grid defined, Lanczos interpolation with
filter size parameter $a = 3$ on the input image is used to determine
the values to assign to the output pixel location.  The output
locations are shifted by 0.5 pixels to let the interpolation select
the value that would be assigned to the center of the output
pixel. This process is repeated for all grid boxes, for all input
images, and for each output image product: the science image, the
variance, and the mask. \czwdraft{The jacobian is multiplied to the
  image value, and squared and multiplied to the variance.  I don't
  understand that.}

As the interpolation constructs the output pixels from more than one
input pixel, there is a covariance term that is must be included.  For
each locally linear grid box, the covariance is calculated from the
kernel in the center of the 128 pixel range.  Once the image has been
fully populated, this set of individual covariance matrices is
averaged to create the final covariance for the full image.

An output catalog is also constructed from the full exposure input
catalog, including only those objects that fall on the warped image.
These detections are transformed to match the new image location, and
to scale the position errors based on the new orientation.

The output image also contains header keywords SRC\_0000, SEC\_0000,
MPX\_0000, and MPY\_0000 that contain the mappings from the warped
pixel space to the input image.  The SRC keyword lists the input OTA
name, and the SEC keyword lists the image section corresponding to the
locally linear grid box.  The MPX and MPY contain the transformation
parameters for the locally linear grid.  \czwdraft{Is this accurate?}

% Read all images and astrometry
% Check which input images overlap with output image. => 8007 when the inputs don't overlap.
% Loop over each image.
% Append detections from input into output detection list.
% Determine transform back from warp pixels to source image.
%% 2nd order polynomial in both x and y for this transformation. and save to header
% Break warp image into 128x128pixel locally linear areas
% Mask non finite pixels as saturated.
% Define Lanczos-3 interpolation over the input image.
% Iterate over warp pixel space (on each locally linear grid) and map interpolated input pixel positions onto warp.
% Warp pixel space is defined as center based, so that's where the intpolation point comes from.
% Covariance calculated based on the interpolation kernel at the center of the ll grid.
% image = interp_image * jacobian
% var   = interp_var * jacobian**2
% mask  = interp_mask
% jacobian = abs(mapXx * mapYy - mapYx * mapXy)
% I don't understand that jacobian.
%


\section{Stacking}

Once individual exposures have been warped onto a common projection
system, they can then be combined without that added concern.  In
order to obtain detections of faint images, and to provide a static
sky image without transient features, we coadd the individual warps
into a stacked image.  Creating this stack also allows a complete
image to be constructed that does not have regions masked due to
falling between devices.

The stacked image is comprised of all warp frames for a given skycell
in a single filter.  The source catalogs and images are loaded into
the \textbf{ppStack} program to do prepare the inputs and stack the
frames while rejecting bad pixels.

Once all files are ingested, the first step is to measure the size and
shapes of the input image PSFs.  We exclude images that have a PSF
FWHM greater than 10 pixels, as those images have the worst seeing and
would degrade the final output stack.  A target PSF for the stack is
constructed from the envelope of all input PSFs, which sets the target
PSF at the largest value among the input PSFs for all radii.  This PSF
is then circularized to prevent any of the input images from being
deconvolved when matched to the target.

The input images also need to be normalized to prevent differences in
seeing and sky transparency from causing discrepancies during pixel
rejection.  From the calibrated input catalogs, we have the
instrumental magnitudes of all sources, along with the airmass, image
exposure time, and zeropoint.  All output stacks are calibrated to a
zeropoint of 25.0 in all filters, and to have an airmass of 1.0.  The
output exposure time is set as the sum of the input exposure times.
With this information, we can determine the relative transparency for
each input image by comparing matched sources between the different
images.  Each image then has a normalization factor defined, equal to
$norm_{image} = (ZP_{image} - ZP_{target}) - transparency_{image} -
2.5 * \log_{10} (t_{target} / t_{image}) - airmassTerm *
(airmass_{image} - airmass_{target})$.  The input source catalog is
adjusted to reflect this normalization, which is also retained for
application when the pixels are combined.  

% PREPARE
% load sources
% load psfs
% determine target as envelope of input psfs
% FWHM clipping at 10
% measure seeing
% -         // M_app = m_inst + zp + c1 * airmass + 2.5log(t) - transparency
%         // EAM : the discussion here was not quite right (or at least sloppy).  Here is a replacement explanation:
%        // For any star, the observed instrumental magnitude on an image and the apparent magnitude are related by:
%        // M_app = m_inst + zp + c1 * airmass + 2.5log(t) - transparency
%        // NOTE the sign of 'transparency'  this must agree with the definition in pmSourceMatch.c. see, eg, line 457 where 
%        // transparency = m_inst + zp + c1 * airmass + 2.5log(t) - M_app 
%        // we want to adjust the input images to be in a consistent flux system so that the
%        // final stack can be generated with a specific target zero point.  Any adjustment to
%        // the flux scale of the image must be made in coordination with the resulting
%        // zeropoint, exposure time, and airmass such that the above relationship yields the
%        // same apparent magnitude for a given star:
%        // m_inst_i : instrumental mags on input image (in)
%        // m_inst_o : instrumental mags on re-normalized image (out)
%        // m_inst_o + zp_o + c1 * airmass_o + 2.5log(t_o) - trans_o = m_inst_i + zp_i + c1 * airmass_i + 2.5log(t_i) - trans_i
%        // m_inst_o = m_inst_i + (zp_i - zp_o) + c1 * (airmass_i - airmass_o) + 2.5log(t_i) - 2.5log(t_o) - trans_i + trans_o
%        // zp_i, airmass_i, t_i, trans_i : reported or measured for input image
%        // zp_o      = zpTarget      (from recipe)
%        // airmass_o = airmassTarget (from recipe)
%        // t_o       = sumExpTime    [sum of input exposure times: once images are scale to this time, they can be avereaged]
%        // trans_o   = 0.0           [obviously!]
%        // we have 2 cases: (a) all reported ZPs are good or (b) some are bad:
%        // (a) FPA.ZP = zp_i + c1 * airmass_i
%        //  --> zp[i] = zp_i + c1 * airmass_i + 2.5log(exptime_i)
%        // (b)  zp[i] = c1 * airmass_i + 2.5log(exptime_i)
%        // NOTE: in case (b), the current code is equating the TARGET zp with the NOMINAL zp, which is wrong.
%        // m_inst_o - m_inst_i = zp[i] - zpTarget - c1 * airmassTarget - 2.5log(sumExpTime) - trans_i

With the normalization factors and target PSF chosen, the convolution
kernels can be calculated for each image.  ISIS kernels are used with
FWHM values of 1.5, 3.0, and 6.0 pixels and polynomial orders of 6, 4,
and 2.  \czwdraft{Skipping this bit because I'm not completely sure I
  understand it.}  The image is then scaled by the normalization as
$renorm = 10^{-0.4 * norm_{image}} / norm_{convolution}$, and the
  variance by the square of that value.  

% MATCH
% match to target PSF.
% use ISIS kernels to do matching/convolution
% Input sources used for psf matching.
% @ISIS.WIDTHS    F32     1.5  3.0  6.0   # Gaussian kernel FWHM values
% @ISIS.ORDERS    S32     6    4    2     # Polynomial orders for ISIS kernels

Once the convolution kernels are defind for each image, they are used
to convolve the image to match the target PSF.  Any input image that
has a $\chi^2$ value larger than 4.0$\sigma$ larger than the median
value is rejected from the stack.

% CONVOLVE
% Normalization to match target zeropoint/exptime
% Reject images with bad match chi^2 values.  MATCH.REJ * rms + median threshold.
% Additional variance from the convolution chi^2
% Calculate image weights based on variance: W_i = 1 / (ROBUST_MEDIAN(variance_image_i) * CovarianceFactor)
% CovarianceFactor = covariance->kernel[0][0]

Following the convolution, and initial stack is constructed.  For a
given pixel coordinate, the values at that coordinate are extracted
from all input images.  Images that have a suspect mask bit (including
the SUSPECT, BURNTOOL, SPIKE, STREAK, STARCORE, and CONV.POOR bit
values) are appended to a suspect pixel list for preferential
exclusion.  Following this, the pixel values are combined and tested
to attempt to identify discrepant values that should be excluded.

If only a single input is available, the initial stack contains the
value from that single input.  If there are only two inputs, the
average of the two is used.  These cases should occur only rarely in
the $3\Pi$ survey, as there are many input exposures that overlap any
particular point on the sky.  The more common case for three or more
inputs constructs a weighted average from the inputs, with the weight
set as a single value for each input image, and defined as the inverse
of the median variance value from that image's associated variance
map.  This weight is used for the image and the exposure weighted
image:

\begin{eqnarray}
  S_{value} &=& \sum_i\left(value_{i} * weight_i\right) / \sum_i\left(weight_i\right) \\
  S_{exp weight} &=& \sum_i \left(exptime_i * weight_i\right) / \sum_i\left(weight_i\right) \\
\end{eqnarray}

The pixel exposure time is simply the sum of the input exposure time values, and the output variance is 

\begin{eqnarray}
  S_{variance} &=& 1 / \sum_i \left( 1 / variance_i \right)
\end{eqnarray}

The output mask value is taken to be zero (no masked bits), unless
there were no valid inputs, in which case the BLANK mask bit is set.

% INITIAL COMBINE
% Calculate weighted mean of input images
% mu = sum_i(f_i * W_i) / sum_i(W_i)
% sigma = 1 / sum_i(1 / W_i)
% nu = sum_i(m_i)
%     // We're not using the input pixel variances to generate a weighted average for the pixel flux (because
%    // that introduces systematic biases), so the variance of the output pixel value should simply be:
%    //     simga^2 = sum(weight_i^2 * sigma_i^2) / (sum(weight_i))^2
%    // This reduces, when the weights are all identically unity, to:
%    //     variance_combination = sum(variance_i) / N^2
%    // and if the variances are all equal:
%    //     variance_combination = variance_individual / N
%    // which makes sense --- the standard deviation of the combination is reduced by a factor of sqrt(N).
% sumValueWeight = sum_i(values * weights)
% sumWeight = sum_i(weights)
% sumVarianceWeight == sum( 1 / variances)
% sumExp  = sum_i(exptimes)
% sumExpWeight = sum_i(exptims * weights)
% mean = sumValueWeight / sumWeight
% var  = 1 / sumVarianceWeight
% exp = sumExp
% expWeight = sumExpWeight

% EXCEPT: if N = 1, accept it.  if N = 2, take average.

%     if (!pmStackCombine(outRO, NULL, stack, maskBad, maskSuspect, maskBlank, kernelSize, iter,
%                        combineRej, combineSys, combineDiscard, useVariance, safe, nminpix, false)) {
%bool pmStackCombine(
%    pmReadout *combined,                // output stacked readout
%    pmReadout *expmaps,                 // output exposure map information
%    psArray *input,                     // input exposures
%    psImageMaskType badMaskBits,        // treat these bits as 'bad'
%    psImageMaskType suspectMaskBits,    // treat these bits as 'suspect'
%    psImageMaskType blankMaskBits,      // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?)
%    int kernelSize,
%    float iter,             0.5
%    float rej,              4.0
%    float sys,              0.1
%    float olympic,          0.2
%    bool useVariance, 
%    bool safe, 
%    int nminpix,
%    bool rejection)
%{

% combineExtract
%% pixels with mask values as suspect are appended to suspect pixel list.
% combinePixels
%% As described above.

Following this initial combination, a ``testing'' loop iterates in an
attempt to identify outlier points.  Again, if only one input is
available, that input is accepted.  If there are two inputs, $A$ and
$B$, then a check is made to see if $(0.5 * (value_A - value_B))**2 >
rej**2 * (variance_A + variance_B + (sys * value_A)**2 + (sys *
value_B)**2)$, where $rej$ is the number of sigma deviant a point
needs to be to be excluded, set to 4.0 for the PV3 processing, and
$sys$ is an estimate of the systematic error, taken to be 0.1.


\czwdraft{This discussion seems out of place, but I'm not sure where a
  better place is.}  Due to the various non-astronomical ghosts that
can occur on GPC1, and the fact that they may not be masked
aggressively enough to ensure all bad pixels are removed, it is
expected that some of the inputs for a given stack pixel are not in
agreement with the others.  In general, there is the population of
input pixel values around the correct astronomical level, as well as
possible populations at lower pixel value (such as due to an
over-subtracted burntool trail) and at higher pixel values (such as
that caused by an incompletely masked optical ghost).  Due to the
observation strategy to image a given field twice to allow for
warp-warp difference images to be constructed to identify transient
detections, higher pixel values that come from sources like optical
ghosts that are a function of pointing will come in pairs as well.
The higher pixel value contaminants are also potentially problematic
as they may appear to be a real source, prompting photometry to be
performed on a false object.  Because of these reasons, there is a
slight preference to reject higher pixel values than lower pixel
values.

If the number of inputs is larger than 6, then a Gaussian mixture
model analysis is run on the inputs to fit two sub populations, and
determine an the likelihood that the distribution is best described by
an uni-modal model.  If this probability is less than 0.05, then the
mean is taken from the bimodal sub population with the largest
fraction of inputs, as this should exclude any sub population
comprised of high pixel value outliers.

If this is not the case (the distribution is likely unimodal) or if
there are insufficient inputs for the mixture model analysis, the
input values are passed to an Olympic weighted mean calculation.  We
set 0.2 as the fraction of the number of inputs to reject through this
process.  This sets the number of bad inputs at $N_{bad} = 0.2 *
N_{input} + 0.5$, where the 0.5 term ensures at least one input is
rejected.  This number is further separated into the number of low
values to exclude $N_{low} = N_{bad} / 2$, which will default to zero
if there are few inputs due to integer arithmatic, and $N_{high} =
N_{input} + N_{low} - N_{bad}$.  After sorting the input values to
determine which values fall into the low and high groups, the
remaining input values have a weighted mean calculated as described
above.

A systematic variance term is necessary to correctly scale how
discrepant points can be from the ensemble mean.  If the mixture model
analysis was run, the Gaussian sigma from the largest sub population
is squared and used.  If this is not available, a 0.1 scaling on the
input values is used.  Each point then has a limit calculated:

\begin{eqnarray}
  limit_{mixture_model} &=& rej**2 * (variance_i + \sigma_{MM}^2) \\
  limit_{default} &=& rej**2 * (variance_i + (0.1 * value_i)**2)
\end{eqnarray}

where $rej$ is the same factor of 4.0 used above.  Each input pixel is
then compared against this limit, and the most discrepant pixel that
has $(value_i - mean)**2$ exceeding this limit is identified.  If
there are suspect pixels in the set those pixels are marked for
rejection, otherwise this worst pixel is marked for rejection.
Following this, the combine and test loop is repeated for a total $0.5
N_{input}$ iterations, or until no more pixels are rejected.

% combineTest
%% if (Ninput > 6) { use KMM }
%% KMM: 
%% Calculate KMMmu KMMsigma KMMpi KMMPunimodal
%% SumWeights = sum(pixelWeights)
%% SysVar = KMMSigma**2 OR (sys * pixelData[i])**2
%% pixelLimts[i] = rej**2 * (pixelVariances[i] + sysVar)
% Iterate 0.5 * Ninput times (at least once)
%% Ninput = 1 => accept
%% Ninput = 2 => if (0.5 * (A - B))**2 > rej**2 * (pixelVariance[A] + pixelVariance[B] + (sys * A)**2 + (sys * B)**2)
%%               then if (suspect) mark reject else mark inspect
%% Else       => if (useKMM and Punimodal < 0.05) median = KMMmean
%%            => else median = combinationWeightedOlympic{}
%%            => if (pixelData - median)**2 > pixelLimits[i] then find single worst deviant pixel value
%% then       => if suspect (mark reject) else (mark reject worst deviant pixel value)


%% combinationWeightedOlympic =>
%% numBad = frac * Ninput + 0.5
%% low = numBad / 2, high = low + numGood - numBad
%% sort(values) => 
%% if (i > low && i <= high) { sumValues = sum_i(values * weights); sumWeight = sum_i(weights)
%% return (sumValues / sumWeight)

% obtain lists of inspect and reject pixels.

% normalize:?
%            float normalise = powf(10.0, -0.4 * norm->data.F32[i]); // Normalisation
%            psBinaryOp(ro->image, ro->image, ``*'', psScalarAlloc(normalise, PS_TYPE_F32));
%            psBinaryOp(ro->variance, ro->variance, ``*'', psScalarAlloc(PS_SQR(normalise), PS_TYPE_F32));

With the initial list of rejected pixels generated, a rejection mask
is made by constructing an empty image that has the rejected pixels
set to a value of 1.0.  This image is then convolved with a 5 pixel
FWHM 0-order ISIS kernel.  Any pixels that are above the threshold of
0.5 are marked as bad and will be rejected in the final convolution.
If more than 10\% of all pixels from an input image are rejected, then
that entire image is rejected as well.

% PIXEL REJECTION
% Construct 15-pixel wide ISIS kernel with 5 pixel FWHM 0-order.
% Construct image of pixels to inspect and convolve with kernel (normalize out kernel power)
% Determine pixels are bad if they're larger than THRESHOLD.MASK = 0.5.
% If more than IMAGE.REJ = 0.1 fraction of pixels are rejected, the entire image is rejected.


\czwdraft{I'm not entirely sure why we do what appears to be a similar
  operation twice.  It also seems odd that this is in the CombineFinal
  step, and not in the Reject step.}  Finally, the rejected pixels are
allowed to grow to include pixels that are neighbors to many rejected
pixels.  The ISIS kernel used in the previous step is used to
determine the largest square box that contains under the limit of
$0.25 * \sum_{x,y} kernel**2$.  This box is then convolved with the
rejected pixel mask to reject their neighbors.

This final list of rejected pixels is passed to the final combination
pass, which does not iterate, and simply excludes the rejected
pixels. \czwdraft{This is a bad paragraph.}

\czwdraft{We make the stacked image, the stacked variance, the stack
  mask, the exposure time mask, the exp weight containing the weighted
  exposure times, and a number image, containing the number of inputs
  used for each pixel.}





% FINAL COMBINE
% Grow rejected pixels 
%% set threshold of (POOR.FRACTION = 0.25) * sum(kernel)**2
%% Choose the largest square box that contains just under that threshold.
%% Convolve that box with the rejected pixels to grow them.
% Run combination pass again, but without doing rejection, simply applying the rejection lists already calculated.
% ::
%      if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, true, true, true)) {
% iter = 0
% combineRej = NAN
% combineSys = NAN
% combineDiscard = NAN
%    if (!pmStackCombine(outRO, expRO, stack, maskBad, maskSuspect, maskBlank, 0, iter, combineRej,
%                        combineSys, combineDiscard, useVariance, safe, nminpix, rejected)) {


The convolved stack products are not retained, as the convolution
reduces the resolution of the final image.  Instead, we apply the
normalizations and rejected pixel maps generated from the convolved
stack process to the original unconvolved input images.  This produces
an unconvolved stack that has the optimum image quality possible from
the input images.  Not convolving does mean that the PSF shape changes
somewhat across the image, as the different FWHM of the input images
print through in the different regions in which they have contributed
to the final image.

% UNCONVOLVED IMAGE
%         if (!ppStackCombineFinal(stack, options->origCovars, options, config, false, true, false, true)) {
% no grow

% only retain unconvolved products.


One benefit of producing the final stacked image from the weighted
mean of the unrejected input images is that faint sources do not have
their contribution removed as much as median filtering would allow.
\czwdraft{I did something to prove this once, but can't find it right
  now.  Comparing the ppStack output catalog to one constructed from a
  simple median filtered stack shows that the ppStack catalog detects
  sources up to 0.XX magnitudes fainter than the median stack.  This
  does increase the possibility of false positives.}


\section{Discussion}

\section{Conclusion}

\end{document}


% http://svn.pan-starrs.ifa.hawaii.edu/trac/ipp/wiki/GPC1_Detrend_Documentation
