Infrared Imaging Data Reduction Software and Techniques
C.N.\ Sabbey
R.G.\ McMahon
J.R.\ Lewis
\& M.J.\ Irwin
% Leave intact We describe the InfraRed Data Reduction (IRDR) software package, a small ANSI C library of fast image processing routines for automated pipeline reduction of infrared (dithered) observations. We developed the software to satisfy certain design requirements not met in existing packages (e.g., full weight map handling) and to optimize the software for large data sets (non-interactive tasks that are CPU and disk efficient). The software includes stand-alone C programs for tasks such as running sky frame subtraction with object masking, image registration and coaddition with weight maps, dither offset measurement using cross-correlation, and object mask dilation. Although we currently use the software to process data taken with CIRSI (a near-IR mosaic imager), the software is modular and concise and should be easy to adapt/reuse for other work. IRDR is available from anonymous ftp to ftp.ast.cam.ac.uk in pub/sabbey.
paper.tex
%ADASS_PROCEEDINGS_FORM%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SAMPLE1.TEX -- ADASS X (2000) ASP Conference Proceedings sample
% paper with minimal markup. Based on the sample from ADASS IX (99).
%
% This is a simple example.  If you want to see a more comprehensive
% sample paper,  take a look at sample2.tex.
%
% Much of the input will be enclosed by braces (i.e., { }).  The
% percent sign, "%", denotes the start of a comment; text after it
% will be ignored by LaTeX.  You might also notice in some of the
% examples below the use of "\ " after a period; this prevents LaTeX
% from interpreting the period as the end of a sentence and putting
% extra space after it.   
% 
% You should check your paper by processing it with LaTeX.  For
% details about how to run LaTeX as well as how to print out the User
% Guide, consult the README file.  (The User Guide is also available
% on-line via 
% http://hea-www.harvard.edu/ADASS
% in the Proceedings section)
%
% If you do not have access to the LaTeX software or a laser printer
% at your site, you can still prepare your paper following the
% instructions in the User Guide.  In such cases, the editors will
% process the file and make any necessary editorial adjustments.
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
\documentstyle[11pt,adassconf]{article}  % Leave intact

\begin{document}   % Leave intact

%-----------------------------------------------------------------------
%			    Paper ID Code
%-----------------------------------------------------------------------
% Enter the proper paper identification code.  The ID code for your
% paper is the session number associated with your presentation as
% published in the official ADASS 2000 conference proceedings.  You can
% find this number locating your abstract in the printed proceedings
% that you received at the meeting or on-line via
% http://hea-www.harvard.edu/ADASS; the ID code
% is the letter-number sequence proceeding the title of your
% presentation.  
%
% This will not appear in your paper; however, it allows different
% papers in the proceedings to cross-reference each other.  Note that
% you should only have one \paperID, and it should not include a
% trailing period.

\paperID{P2-43}

%-----------------------------------------------------------------------
%		            Paper Title 
%-----------------------------------------------------------------------
% Enter the title of the paper.

\title{Infrared Imaging Data Reduction Software and Techniques}

%-----------------------------------------------------------------------
%		          Authors of Paper
%-----------------------------------------------------------------------
% Enter the authors followed by their affiliations.  The \author and
% \affil commands may appear multiple times as necessary.  List each
% author by giving the first name or initials first followed by the
% last name.  Authors with the same affiliations should grouped
% together. 
%
% Try to limit the front matter to no more than three \author
% commands.  Group authors with the same affiliations.  Too many
% \author commands fills the first page of the paper with little
% actual text.

\author{C.N.\ Sabbey, R.G.\ McMahon, J.R.\ Lewis, \& M.J.\ Irwin}
\affil{Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK}

%-----------------------------------------------------------------------
%			 Contact Information
%-----------------------------------------------------------------------
% This information will not appear in the paper but will be used by
% the editors in case you need to be contacted concerning your
% submission.  Enter your name as the contact along with your email
% address.

\contact{Chris Sabbey}
\email{sabbey@ast.cam.ac.uk}

%-----------------------------------------------------------------------
%		      Author Index Specification
%-----------------------------------------------------------------------
% Specify how each author name should appear in the author index.  The 
% \paindex{ } should be used to indicate the primary author, and the
% \aindex for all other co-authors.  You MUST use the following
% syntax: 
%
% SYNTAX:  \aindex{LASTNAME, F. M.}
% 
% where F is the first initial and M is the second initial (if
% used).  This guarantees that authors that appear in multiple papers
% will appear only once in the author index.  

\paindex{Sabbey, C. N.}
\aindex{Lewis, J. R.}     
\aindex{McMahon, R. G.}     
\aindex{Irwin, M. J.}

%-----------------------------------------------------------------------
%			Subject Index keywords
%-----------------------------------------------------------------------
% Enter up to 6 keywords describing your paper.  These will NOT be
% printed as part of your paper; however, they will be used to
% generate the subject index for the proceedings.  There is no
% standard list; however, you can consult the indices for past ADASS
% proceedings (http://iraf.noao.edu/ADASS/adass.html). 

\keywords{infrared data reduction, image processing}

%-----------------------------------------------------------------------
%			       Abstract
%-----------------------------------------------------------------------
% Type abstract in the space below.  Consult the User Guide 
% (http://hea-www.harvard.edu/ADASS)
% for a list of supported macros (e.g. for typesetting special
% symbols). 

\begin{abstract}          % Leave intact

We describe the InfraRed Data Reduction (IRDR) software package, a small
ANSI C library of fast image processing routines for automated pipeline
reduction of infrared (dithered) observations.  We developed the software
to satisfy certain design requirements not met in existing packages (e.g.,
full weight map handling) and to optimize the software for large data sets
(non-interactive tasks that are CPU and disk efficient).  The software
includes stand-alone C programs for tasks such as running sky frame
subtraction with object masking, image registration and coaddition
with weight maps, dither offset measurement using cross-correlation,
and object mask dilation.  Although we currently use the software to
process data taken with CIRSI (a near-IR mosaic imager), the software
is modular and concise and should be easy to adapt/reuse for other work.
IRDR is available from anonymous ftp to ftp.ast.cam.ac.uk in pub/sabbey.

\end{abstract}

%-----------------------------------------------------------------------
%			      Main Body
%-----------------------------------------------------------------------
% Place the text for the main body of the paper here.  You should use
% the \section command to label the various sections; use of
% \subsection is optional.  Significant words in section titles should
% be capitalized.  Sections and subsections will be numbered
% automatically. 

\section{Introduction}

The Cambridge Infrared Survey Instrument (CIRSI) is a near-IR mosaic
imager that contains a 2 x 2 array of Rockwell Hawaii I 1024 x 1024
detectors (Beckett et al.\ 1996; Mackay et al.\ 2000).  CIRSI has been in
operation for about two years and we have obtained almost 1 Terabyte of
imaging data.  The uniquely wide field accessible by CIRSI on a 2--4m
class telescope makes it ideal for large-area surveys to moderate depth.
Two such surveys are currently in progress with preliminary results
including the measurement of galaxy clustering at intermediate redshift
(McCarthy et al.\ 2001), and the demonstration of a reddening-independent
quasar selection technique based on combined deep optical and near-IR
color diagrams (Sabbey et al.\ 2001).

However, the CIRSI data reduction poses several challenges.  With the
large data rate (5--10 GB of data taken per night, $\sim 100$ nights per
year, and a significant data backlog currently) the software has to be
very efficient and completely automated.  Also, the software should
handle diverse data sets, from Galactic center observations to very
sparse fields at high Galactic latitude.  We generate wide-field, deep
mosaic images from thousands of individual images taken over many nights.
With the gaps between detectors comparable to detector size, we make
filled mosaic images using the coadded dither sets from different chips
and telescope pointings.  Thus weight maps and accurate astrometry are
crucial and the common simplification of clipping the dither sets to
their intersection region is not appropriate.

Although we use existing software packages when possible (see below),
we decided to write the core image processing routines ourselves to
satisfy certain design requirements.  For example, we wanted a two pass
reduction with object masks derived from the first pass coadded dither
sets, subpixel image registration and coaddition that uses full weight
maps and does not clip images, optimizations for CPU/disk efficiency,
customized artifact cleaning (destriping and defringing), and reusable
tools (from having stand-alone tasks to glue together using a high
level scripting language, to making C library calls, or even extracting
portions of source code).  The basic processing steps, described below,
are: flatfield correction, running sky frame subtraction, dither offsets
measurement, dither set coaddition, and mosaic image creation.

% as stand-alone tasks to run from some scripting language, to making C 
% library calls, or even extracting portions of source code
%at the source code level or for gluing together tasks with some high
%level scripting language).  
% Not a general purpose image processing suite like Eclipse or IRAF, but
% rather concentrates on doing certain specific tasks well.

\section{Data Reduction}
\subsection{Flatfield Correction}
Flatfield images are produced by subtracting the stack median of lamp-off
domeflats from the stack median of lamp-on domeflats.  The flatfield
images are divided by the mode of the chip 1 flatfield to produce a gain
map per chip.  Bad pixels are automatically identified in the gain maps
(and set to 0.0) by looking for outliers ($> 5\sigma$ from the median
in 15x15-pixel blocks) or pixels with extremely low or high sensitivity
($> 30$\% from the median gain).  Because bad pixels often occur in
clumps we eliminate bad pixels during image coaddition rather than by
interpolation in an initial cleaning pass.  The data frames are multiplied
by the inverse of the gain map.

The image stacking is done using {\tt cubemean.c}, which calculates the
median, robust standard deviation, or robust mean plane with a choice
of weights (none, scalar, or maps) and image scaling or zero offsets
using the image modes.  This is not a general purpose tool like IRAF's
{\tt imcombine}, but for the specific task of calculating the median plane
was found to be 2.5 times faster (for a stack of 50 of our data frames).
The flatfield image is converted to a gain map and bad pixels identified
using {\tt gainmap.c}.  The data are flatfield corrected using {\tt flat.c}.

\subsection{Running Sky Frame Subtraction}
For each data frame, a sky image is constructed from the robust mean
of the 8 nearest frames in the observation sequence and subtracted
({\tt skyfilter.c}).  Objects detected in the coadded dither sets from
the first pass reduction are masked out during sky frame creation in the
second pass.  The object masks are produced using the checkimage {\tt
OBJECTS} option to SExtractor (Bertin and Arnouts 1996), which produces
a FITS image with non-object pixels set to 0.  This is simpler and more
effective than building masks from a catalog of object coordinates and
shape parameters.  The object regions (detection isophotes) generated by
SExtractor are then expanded by a multiplicative factor of 1.5 (using
{\tt dilate.c}), thereby growing the mask regions for large objects
more than small objects.  The object masks for individual frames are
obtained on the fly using pixel offsets (i.e., the dither offsets)
into the master dither set masks.

% (and merged with weight maps also generated on the fly).
% , and not used in calculating the sky image.  

% also, non-coadded images are stored and read/written as short integers
% (2 bytes deep) instead of floating point (4 bytes)
% advantage of input short integer data to do fast statistics calculations.
% ten times faster than optimized median
% functions for our data images even faster than byte swapping.

Running sky frame subtraction is normally a significant bottleneck in
processing infrared imaging, so optimizations are important.  To do
running sky frame subtraction for a stack of N images, the program {\tt
skyfilter.c} uses a sliding window (circular buffer of image pointers) to
only require N image reads and N image mode calculations.  In contrast,
putting this logic into a script normally involves N $\times$ M image
reads and mode calculations, where M is the width of the sky filter
in frames.  Also, the disk I/O (and storage) for non-coadded data uses
short integers (2 bytes deep), even though most calculations are done
in floating point (4 bytes).  Some calculations work with short integer
data however to allow optimizations.  For example, the almost trivial
distribution sort can be used to obtain an image histogram, sorted
image array, and median value in $O(n)$ time ($\approx 5$ times faster
than running an optimized median routine on our data images).

%\begin{verbatim}
%for (i = 0; i < n; i++)  /* fill histogram using image data */
%    hist[data[i]]++;
%
%for (j = 0, i = 0; j < m; j++)     
%    while(hist[j]--)   /* unstack histogram into sorted array */
%        data[i++] = j;
%\end{verbatim}
%
%If the cumulative histogram and cumulative histogram weighted by bin
%number are then calculated, 
%
%Then binary searches into the cumulative histogram can be used to do
%iterative mode and sigma calculations.  faster than byte swapping.
%
%Object mask dilation is done using dilate.c
%With C program we only need N reads and N mode calculations for stack of N 
%images, whereas IRAF scripts do N * M where M is the number of frames in
%the sky filter.  Sorting/Median/Mode calculations can be done faster than
%byte swapping if we use the fact that the data are short integers 
%(distribution sort etc).

\subsection{Dither Offsets Measurement}
A typical dither sequence consists of nine observations in a $3 \times
3$ grid with offsets in each direction of $\approx 10$ arcsec.  The
approximate dither offsets stored in the FITS header WCS information are
refined using cross-correlation analysis ({\tt offsets.c}).  The non-zero
(object) pixels of the reference frame object mask (SExtractor OBJECTS
image) are stored in a pixel list (x, y, brightness), and this list is
cross-correlated against the object mask images of the following frames
in the dither set.  The SExtractor OBJECTS image conveniently removes
the background (important for cross-correlation methods) and identifies
the object pixels more reliably than a simple thresholding algorithm
(e.g., especially in images with a non-flat background, large noise, and
cosmic rays).  Using an object list in the cross-correlation focuses on
the pixels that contribute to the cross-correlation signal and is faster
than cross-correlating two images.

The cross-correlation technique uses coordinate, magnitude, and shape
information, and was found to be more reliable than matching object
coordinate lists (the improvement was noticed in extreme cases, like
Galactic center images and nearly empty fields with an extended galaxy).
A subpixel offset measurement accuracy of about $0.1$ pixels is obtained
by fitting a parabola to the peak of the cross-correlation image.
In terms of speed, this cross-correlation method was found to be $\approx
10$ times faster (for typical survey data and a relatively large search
box of 100 pixels) than IRAF STSDAS {\tt crosscor}.  Although the success
rate is $\approx 100$\%, failure is indicated by an offset measurement
corresponding exactly to the border of the search area, or a small fraction
of object pixels overlapping in the aligned data images.

\subsection{Dither Set Coaddition}
A weight map is generated on the fly for each data frame with the weight
for pixel $p_{i}$ given by: $w_{i} = g_{i} \cdot t / V$, where $g_{i}$ is the
gain for pixel $p_{i}$ (0.0 for bad pixels), $t$ is the exposure time, and $V$
is the image variance.  The data frames and corresponding weight maps in
the dither set are registered using bi-linear interpolation modified to
account for bad pixels and image weights.  Each output (interpolated) pixel 
value $P$ is calculated from the weighted average of the four overlapping
input pixels $p_{i}$ of the input image:
$$
P = \frac{1}{W} \sum_{i = 1}^{4} a_{i} w_{i} p_{i}, \; \; \; \;
W = \sum_{i = 1}^{4} a_{i} w_{i}
$$
where $a_{i}$ are the pre-calculated fractional areas of overlap of $P$
with $p_{i}$, and $w_{i}$ are the image weights for $p_{i}$.  The weight
maps are registered similarly but using a weighted sum to calculate
the weight $W$ for pixel $P$.  A different registration method that
people sometimes recommended is to replicate each pixel of an image
into N x N pixels and do an integer shift in units of these new pixels.
Although this will approximate the above bi-linear interpolation scheme
as N becomes large, it requires more work for less precision.  We have
not tested higher order interpolators (e.g., Devillard 2000), although
our default approach is fast and reasonable, especially given the low
signal-to-noise ratio of our individual data frames.

%fast and reasonable, especially for our individual background-limited,
%noisy images.

%especially considering the low signal-to-noise ratio of our individual
%data frames.
%reasonable, especially for noisy images.

The dither frames are combined by calculating the weighted mean pixel
value at each $(x,y)$ position of the dither stack, with pixel values $>
5 \sigma$ from the median at each position rejected.  Images borders are
added during registration to avoid clipping the data to the intersection
of the dither frames.  The standard deviation ($\sigma$) at each position
is calculated from $\sigma = {\rm MAD} / 0.6745$, where MAD is the median
absolute deviation from the median (we do not use a simpler method such
as minmax rejection because we often take averages of small numbers of
values, e.g., 5--9 frames per dither set).  The weight maps are combined by
calculating the sum at each $(x,y)$ position of the stack of weight maps
(for pixels not clipped during coaddition).  The program used for the
above is {\tt dithercubemean.c}.

\subsection{Mosaic Creation}

The current astrometry pipeline is a small Perl script that
runs SExtractor to produce an object catalog for each coadded
dither set, runs APMCAT (a stand-alone C program available from
www.ast.cam.ac.uk/$\sim$apmcat/) to download over the network the APM
sky coordinates of objects in each field of view, then runs IMWCS from
WCSTools (Mink 1999) to calculate the astrometry fit and update the WCS
information in the FITS headers.  The coadded dither sets and weight
maps from different chips, telescope pointings, and nights are then
drizzled onto a wide-field mosaic image using EIS Drizzle (available
from www.eso.org/eis).  The astrometry residuals between the mosaic
image and the APM catalogue show a random error of $\sigma \approx 0.3$
arcsec without significant systematic effects.

%\acknowledgments

%-----------------------------------------------------------------------
%			      References
%-----------------------------------------------------------------------
% List your references below within the reference environment
% (i.e. between the \begin{references} and \end{references} tags).
% Each new reference should begin with a \reference command which sets
% up the proper indentation.  Observe the following order when listing
% bibliographical information for each reference:  author name(s),
% publication year, journal name, volume, and page number for
% articles.  Note that many journal names are available as macros; see
% the User Guide
% (http://hea-www.harvard.edu/ADASS) for a listing "macro-ized" journals.   
%
% Note the following are some of the tricks that can be used:
%
%   o  \& is used to format an ampersand symbol (&).
%   o  \'e puts an accent grave over the letter e.  See the User Guide
%      for details on formatting special characters.  
%   o  "\ " after a period prevents LaTeX from interpreting the period 
%      as an end of a sentence.
%   o  \aj is a macro that expands to "Astron. J."  See the User Guide
%      for a full list of journal macros
%   o  \adassvii is a macro that expands to the full title, editor,
%      and publishing information for the ADASS VII conference
%      proceedings.  Such macros are defined for ADASS conferences I
%      through IX.
%   o  When referencing a paper in the current volume, use the
%      \adassix and \paperref macros.  The argument to \paperref is
%      the paper ID code for the paper you are referencing.  See the 
%      note in the "Paper ID Code" section above for details on how to 
%      determine the paper ID code for the paper you reference.  
%
\begin{references}
\reference Beckett, M.G., et al.\, 1996, SPIE, 2871, 1152
\reference Bertin, E., and Arnouts, S. 1996, A\&AS, 117, 393
\reference Devillard, N. 2000, The Messenger, 100, 48
\reference Mackay, C.D., et al.\ 2001, SPIE, 4008, 1317
\reference McCarthy, P., et al.\ 2001, in ESO ``Deep Field'' Workshop, held in Garching, Germany, October 2000, astro-ph/0011499
\reference Mink, D., 1999, in ASP Conf. Ser. 172, ADASS VIII, eds. D.M. Mehringer, R.L. Plante \& D.A. Roberts (San Francisco: ASP), 498
\reference Sabbey, C.N., et al. 2001, in New Era of Wide Field Astronomy, held in Preston, England, August 2000, astro-ph/0012294
\end{references}

% Do not place any material after the references section

\end{document}  % Leave intact
