%%% $Id: imageCombination.tex,v 1.2 2005-04-08 19:51:57 price Exp $
\documentclass[panstarrs]{panstarrs}

% basic document variables
\title{Image Combination Algorithm}
\subtitle{A Recommendation}
\shorttitle{Modules SDRS}
\author{Paul Price}
\audience{Pan-STARRS PMO}
\group{Pan-STARRS Algorithm Group}
\project{Pan-STARRS Image Processing Pipeline}
\organization{Institute for Astronomy}
\version{DR}
\docnumber{PSDC-430-???}

\setlength{\topsep}{-2pt}
  
\begin{document}
\maketitle
\sloppy

% -- Revision History --
% provide explicit values for the old versions
% use '\theversion' for the current version (set above)
% use \hline between each table row
\RevisionsStart
% version  Date            Description
DR & 2004 Dec 09 & Draft \\ \hline
\RevisionsEnd

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

\DocumentsInternal
PSDC-430-011  &   Pan-STARRS PS-1 IPP System/Subsystem Design Description \\ \hline
\DocumentsExternal
Posix Standard & Open Group Based Specifications Issue 6, IEEE Std 1003.1, 2003 \\
\DocumentsEnd

\tableofcontents
\pagebreak 
\pagenumbering{arabic}

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

\section{Introduction}

The goals of the \PS{} system can be (overly?) simply expressed
as the following:
\begin{enumerate}
\item To make deep images of the static sky; and
\item To identify variable/moving objects
\end{enumerate}
For the prototype telescope, PS1, multiple images will be taken
consecutively, while for the full \PS{} system, PS4, multiple
images will be taken simultaneously.  In either case, there is a need
to combine multiple images.

Cosmic-ray (CR) hits on the detectors have the potential to influence
both of the goals of the \PS{} system through corrupting certain
pixels on the detectors.  The addition of a pixel influenced by a CR
hit into the static sky image can result in apparent objects not
really present on the sky, and affect the photometry of real objects.
Consequently, CRs will also produce false detections of apparently
variable objects.  Hence it is necessary to flag pixels influenced by
CRs so that they are not propagated to the combined image from the
input images.

The purpose of this document is to propose an algorithm suitable for
combining multiple images in the \PS{} system.

\section{Requirements}

We assume at the outset that multiple images (motivated by the full
\PS{} system concept, we invariably adopt four input images in this
study) have been obtained and photometrically calibrated such that
counts are proportional to flux from the sky, with the relative
scalings between the images known, along with estimates of the Poisson
error for each pixel (which may be calculated from the scaling, gain,
readnoise and background).  In addition, we also assume that the
images have been astrometrically calibrated, such that a
transformation for each image to (and from) a common output coordinate
system is known.  This is the case in the \PS{} Image Processing
Pipeline (IPP) at the commencement of the Phase 4 processing.  It is
worth mentioning that we expect that the input images have similar
seeing, since no effort is budgeted in this step for PSF matching.

The task before us is to take these images and combine them (which, in
addition to simply combining, also involves transforming the images to
the common output coordinate system), rejecting those pixels from the
combination which are affected by CRs, while leaving untouched those
pixels which remain unaffected by CRs.  In particular those pixels
which comprise the cores of stars (even stars which are, to some
degree, undersampled) should not be rejected.

What are the specific requirements for this image combination?  Of
course, there will always be a trade-off between the fraction of CRs
on the image which will be identified and rejected, and the fraction
of ``clean'' pixels on the image which are falsely identified as CRs
--- such a trade-off is always unavoidable since there is noise in the
input images, and since the input images are populated by stars which
may be (as will be the case for \PS{} in excellent seeing)
undersampled.  We propose the following requirements:
\begin{enumerate}
\item {\bf Given four input images with stars at a density of 400 per
  512$^2$ pixels with total flux distributed as $N(m<M) \propto
  10^{0.35M}$ down to 21.5 mag (using a magnitude zero point of
  25~mag) with a FWHM of 2~pixels, the fraction of detector pixels
  that are falsely rejected as CRs shall not exceed 1\%.}  This value
  corresponds to an order of magnitude less than an estimate of the
  fraction of the focal plane not covered by working detectors, so
  that the rejection of pixels in the image combination will not be
  the major determining factor in the working detector area obtained
  in each exposure.  The seeing here is the median seeing expected for
  PS4 (0.6'' with 0.3'' pixels), the limiting magnitude and magnitude
  zero point are appropriate for PS1, the stellar density is that
  expected for $b \sim 10^\circ$ down to $V=22$~mag ($2\times10^5$ per
  square degree) which we obtain from the PS1 DRM, and the apparent
  magnitude distribution corresponds to that expected in the AP
  survey.

\item {\bf Given four input images which overlap in the output
  coordinate frame (having the same scale as the input coordinate
  frames, and random, subpixel shifts), and are completely devoid of
  stars, consisting only of a flat background with Gaussian noise of
  standard deviation $\sigma$ and cosmic rays (or hot pixels) with
  (single pixel) flux (added to the background) equally distributed
  from 0 to $20\sigma$ covering a fraction of 0.7\% of each input
  image, the fraction of pixels in the output image with a flux
  greater than 1.5$\sigma$ above the background level shall not exceed
  0.2\%.}  This follows a Gaussian distribution, since the standard
  deviation expected in the combined image would be $\sigma/2$, and
  hence we expect 99.8\% of pixels to lie within $3\sigma/2$ if the
  CRs have been properly rejected.  The Gaussian distribution serves
  as an approximation to Poisson noise from the sky.  The number of CR
  pixels are chosen to match the number of CR pixels measured on an
  OPTIC dark (0.62\%; this value should be diluted when sky noise is
  added), while the CR flux distribution is chosen to be biased
  towards low flux events which are more difficult to detect than
  bright events; consequently, this is a strong test of the
  algorithm's ability to identify and mask CRs.
\end{enumerate}

To these, we may add the requirement that the implementation of any
proposed algorithm must satisfy the throughput requirements of the
IPP.


\section{Approach}

CRs (along with other artifacts that we might seek to reject) have two
important properties that we may use to identify them in images:
\begin{enumerate}
\item They have sharp boundaries instead of a stellar PSF.
\item They are (generally) not present at the same location on
multiple images of the same piece of sky.
\end{enumerate}
Using both of these methods together should result in clean combined
images.  Relying solely on measuring the gradients in order to
identify sharp boundaries would generate difficulties in combining
critically sampled and under sampled images.

One desirable feature would be to mask CR pixels on the detector
before transforming the image.  The reason for this is that mapping
from the detector coordinate system to the output coordinate system
smears out cosmic rays, so that an event that affected only a single
pixel on the detector may affect multiple pixels on the output image.
The result of this is that CRs in the output image consist of a core
of high flux plus wings in the surrounding pixels.  While the core may
be easily identified through comparison with other images, the wings
are harder to pick up.  Some combination methods (e.g., Tonry's {\tt
autoclean}) attempt to identify the wings by using a lower threshold
in the combination for pixels near pixels already identified as CRs,
but a cleaner solution would be to mask the CR before mapping.
Consequently, we propose an algorithm that transforms the input
images, combines with rejection, identifies those pixels on the source
images that correspond to the rejected pixels in the combination, and
mask those before repeating the combination.  This process is used to
reduce images from HST using Drizzle/Multidrizzle.

The algorithm proposed here, therefore, maps the input images to the
output coordinate system, and combines the transformed images with
rejection.  Then, for each of the input images, a map is made of those
pixels that were rejected in the combination (set to one if the pixel
was masked, zero otherwise), and this map is transformed back to the
coordinate system of the original input image.  Pixels in this
transformed map that exceed a specified threshold (\code{frac}), and
for which the local gradient is less than another specified threshold
(\code{grad}) are masked as cosmic rays.  The impact of these maskings
on the combined image is then propagated through the transformation
and combination (with no rejection) stages to yield a clean, combined
image.

A short comment about the calculation of the local gradient is
necessary.  Simply calculating the gradient from the image in which a
pixel is suspected of being a CR can be problematic.  Firstly, CRs are
not always single pixel events, but often consist of streaks in which
several pixels are affected, which can bias the calculation of the
gradient.  And secondly, in the case where the seeing yields
critically sampled (or even under sampled) images, the gradient in the
core of a star can be quite high, causing it to be falsely identified
as a CR.  We avoid these problems by calculating the median flux of
neighboring pixels on the other input images (using the known
coordinate transformations), and using the mean difference of these
from the flux of the suspect CR pixel.  Because we look at the other
images, which should be devoid of CRs at the corresponding position,
the first problem is circumvented; and because of the slight
(sub-pixel) shifts between the input images, the gradient can be
downplayed in the case that the pixel is due to the core of a star,
and so the second problem is minimized.

\tbd{Discussion about error images}

\section{Implementation}

We designate our implementation of the above algorithm ``STAC'', for
Simultaneous Telescope Array Combination.  It is built on top of the
\PS{} Library (currently \code{psLib-rel5alpha1}), making use of its
vectors, arrays, images, polynomial transformations, statistics and
tracing functionalities.

The current version does not meet the speed requirements of the IPP,
but we expect that optimising the code for speed should be fairly
simple (currently, the forward--backward--forward methodology
described above is implemented for every pixel of each input image,
but it need only be implemented for pixels which trigger certain
conditions).  We set this requirement aside for the time being,
confident of meeting it in the future, and concern ourselves with the
performance of the algorithm in detecting CRs.

Given input images with no CRs (according to the first requirement
above), STAC falsely identifies $22+109+161+47 = 339$ pixels as CRs in
the four individual $512\times512$ input images, or 0.03\% of the
total number of input pixels.  This meets the requirement by over an
order of magnitude.

Given input images with 2000~CRs (according to the second requirement
above), STAC identifies $1692+1222+1256+1624 = 5794$ pixels as CRs in
the four individual $512\times512$ input images, leaving 29 pixels in
the combined image exceeding 1.5~times the estimated standard
deviation of the input images, or 0.01\% of the pixels in the output
image.  This meets the requirement by an order of magnitude.

Since it meets both of the proposed requirements, STAC (modulo the
speed issue, which should be easily solved) is suitable for use in the
\PS{} IPP.

\section{Alternatives}

Here we examine two alternatives to STAC, as a check to our
recommendation.  The first alternative is a completely naive
combination of the shifted input images using IRAF's \code{imcombine}
(with appropriate rejection).  The second is the more sophisticated
\code{autoclean} developed by John Tonry, which also acts on shifted
input images, but uses a number of rejection thresholds to identify
and reject faint CRs.  We will apply the same tests to these
alternatives as for STAC.

\subsection{imcombine}

We combine the shifted input images using \code{avsigclip} rejection:

\begin{verbatim}
Dec 14 15:19: IMCOMBINE
  combine = average, scale = none, zero = none, weight = none
  reject = avsigclip, mclip = yes, nkeep = 1
  lsigma = 3., hsigma = 3.
  blank = 0.
                Images 
    test_0.fits.shift.1
    test_1.fits.shift.1
    test_2.fits.shift.1
    test_3.fits.shift.1
\end{verbatim}

Given input images with no CRs, \code{imcombine} identifies 5552
pixels as CRs, or 0.5\% of the total number of input pixels.  This
satisfies the requirement by a factor of two.  The combined image
appears noisier than the output image produced by STAC (standard
deviations of 3.6 compared to 3.2).

Given input images with 2000~CRs each, \code{imcombine} leaves 495
pixels in the combined image exceeding the threshold, or 0.19\%.  This
barely satisfies the requirement.  The histogram of pixel values shows
a distinct tail at the high end.

We conclude that STAC performs better in both departments.

\subsection{autoclean}

We combine the shifted input images using the following
\code{autoclean} command line:

\begin{verbatim}
autoclean test.000 4 manual scale=unity zero=zero float average mask trigger=0,2.7,2.0 grad=3 domino=3 eadu=1
\end{verbatim}

Using this setup, we performed the same tests as for the other
combination methods.  In the test on a stellar field with no CRs,
\code{autoclean} falsely identifies 81 pixels as CRs, or an amazingly
low 0.007\%.  In the test on a blank field with CRs added,
\code{autoclean} masks 11286 pixels as CRs, leaving 412 pixels above
the background, or 0.16\% of the area.  As in the case of
\code{imcombine}, this satisfies the requirement, but is still a
factor of 4 above that of STAC; and the histogram still displays a
tail.

Note that the number of pixels masked here cannot be compared directly
with that of STAC, since the latter masks in the source coordinate
fram, while the former masks in the output frame where each CR pixel
from the source occupies multiple pixels.  Dividing by 4 should yield
a more direct comparison, since sub-pixel shifting in both $x$ and $y$
results in (single pixel) CRs being smeared over 4 pixels.  Doing
this, \code{autoclean} masks 2822 CRs in the source, while STAC masks
5875.

We tried other parameters in an attempt to increase the number of
rejected CRs.  Brian Schmidt uses the following as his choice of
default:
\begin{verbatim}
autoclean test.000 4 manual scale=unity zero=zero float average mask eadu=1 scour=50,4.5,2.5
\end{verbatim}
This results in 9862 pixels masked, and 545 pixels above the
threshold.

Combining the \code{trigger} and \code{scour} options:
\begin{verbatim}
autoclean test.000 4 manual scale=unity zero=zero float average mask eadu=1 trigger=0,2.7,2.0 scour=50,4.5,2.5
\end{verbatim}
This results in 11713 pixels masked, and 239 pixels above the
threshold.

Lowering the thresholds for masking pixels:
\begin{verbatim}
autoclean test.000 4 manual scale=unity zero=zero float average mask eadu=1 trigger=0,2.0,1.5 scour=0,1.5,1.0
\end{verbatim}
This results in 12935 pixels masked, and 293 pixels above the
threshold.

The best result we achieved was:
\begin{verbatim}
autoclean test.000 4 manual scale=unity zero=zero float average mask eadu=1 trigger=200,1.5,1.0 scour=200,1.0,0.5
\end{verbatim}
For this, 15419 pixels were masked, resulting in 123 pixels above the
threshold.  Using the same parameters, no pixels are identified as CRs
when using the stellar field devoid of CRs, but this is likely due to
it having a different background level.  Raising the background to the
same level as for the CR image, \code{autoclean} falsely identifies 95
pixels as CRs, whereas STAC identifies 15 for the same field.

We conclude that STAC has slightly better performance than
\code{autoclean}, but \code{autoclean}'s performance is still very
good.

The difference in speed between the two codes is, of course
impressive, with \code{autoclean} being much faster than STAC.  Of
course, STAC's speed is yet to be optimised, and \code{autoclean}
doesn't do the transformations, so the comparison isn't fair.

\section{Conclusion}

STAC and \code{autoclean} both perform well.  STAC performs slightly
better; \code{autoclean} is much faster (by an order of magnitude ---
0.35 seconds versus 3.5 seconds for STAC).  However, there exist
simple optimisations that can be made within STAC which will vastly
improve the speed.  If STAC can be optimised to meet the time budget
for Phase 4, and close to \code{autoclean}, then it should be accepted
on the basis of its superior performance.  If the time budget is
tight, then \code{autoclean} could be used profitably.

\appendix

\section{Addendum: Timing}

Optimisation of the STAC code results in 60 sec to combine a 4k square
image consisting of stars and CRs on \code{alala}, a dual Opteron
2.2~GHz machine.  This appears to be within the time budget for Phase
4.  The initial transformation accounts for over half of this time
(about 37 sec).  \code{autoclean} takes approximately 20 sec
(including I/O) to do the combination.  We conclude that timing is no
longer an issue.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{plain} \bibliography{panstarrs}

\end{document}

