Index: /trunk/doc/misc/imageCombination.tex
===================================================================
--- /trunk/doc/misc/imageCombination.tex	(revision 2794)
+++ /trunk/doc/misc/imageCombination.tex	(revision 2794)
@@ -0,0 +1,389 @@
+%%% $Id: imageCombination.tex,v 1.1 2004-12-23 01:56:58 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 exceeds 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} Release 3), 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 $580+479+57+483 = 1599$ pixels as CRs
+in the four individual $512\times512$ input images, or 0.15\% of the
+total number of input pixels.  This meets the requirement by almost an
+order of magnitude.  As an aside, the low number of false CRs
+identified in the third image appears to be due to the nearly integral
+pixel shift to the output coordinate system.
+
+Given input images with 2000~CRs (according to the second requirement
+above), STAC identifies $1580+1583+1273+1439 = 5875$ pixels as CRs in
+the four individual $512\times512$ input images, leaving 96 pixels in
+the combined image exceeding 1.5~times the estimated standard
+deviation of the input images, or 0.04\% of the pixels in the output
+image.  This meets the requirement by a factor of 5.
+
+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}
+