HST/NICMOS Pedestal Estimation Software
**********************************************************************
* HST/NICMOS Pedestal Estimation & Quadrant Equalization Software *
**********************************************************************
***************************
PART I: GENERAL INFORMATION
***************************
Author
------
* Roeland P. van der Marel,
address : Space Telescope Science Institute
3700 San Martin Drive
Baltimore, MD 21218
Tel : (+1) 410 338 4931
Fax : (+1) 410 338 4596
e-mail : marel@stsci.edu
homepage : http://sol.stsci.edu/~marel/
History
-------
* 1998 : algorithm and program development,
with skilled help and testing by David Zurek.
See Section VI below for details about the various versions of this
software that have been released.
Summary of the method
---------------------
The software reads a calibrated (i.e., flat-fielded) NICMOS image and
also the image that was used for the flat-fielding. It then estimates
the size of any constant bias level (the `pedestal') that may be
present in the data, but was not subtracted before flat-fielding. On
output a new frame is written that has the best-fitting constant
pedestal subtracted. An algorithm that attempts to remove DC
differences between quadrants is also included. If requested, the
program can be used in a mode that only applies the quadrant
equalization algorithm, and omits the pedestal determination.
References
----------
N/A
Acknowledgments
---------------
If you have found this software useful for your research, I would
appreciate an acknowledgment with reference to `use of the Pedestal
Estimation and Quadrant Equalization Software developed by Roeland
P. van der Marel'.
Bug Reports
-----------
Please send bug reports and important comments or questions to
marel@stsci.edu
Caveats
-------
This software comes without guarantee and on an `as is' basis. It is
not being actively maintained or upgraded (however, I will try to
repair bugs if you tell me about them).
I developed this software because I am a GO with NICMOS data. I am not
a member of the NICMOS team at STScI, and this is not officially
endorsed/supported STScI software.
All testing of the code was done on a Sun Ultra 170 workstation with
the Solaris operating system in early 1998.
It is not guaranteed that the software will work without problems on
other platforms, with different operating systems and/or different
compilers. However, the code is close to standard fortran, so if any,
only minor revisions will be necessary.
The programs read and write files in IRAF (.imh) format, that must be
in real format (see Section III below for details). Use of the IRAF
format requires linking of the code to standard IRAF libraries.
This software is freely available. However, I do not encourage people
to distribute this software to third parties. Instead, tell people to
get the latest version of this software from the software section of
my homepage (http://sol.stsci.edu/~marel/software.html).
Accurate determination of a constant bias from an arbitrary image is
difficult. This software provides a possible solution, yet it may not
work as expected in all possible situations. Therefore, please give
proper care to the use of this software, and be sure to carefully
examine both your data and the results of this software. Using this
software as a black-box may lead to erroneous results.
Another point to keep in mind is that this program will, by
definition, search for a constant bias (although this bias may be
different in each quadrant of the NICMOS detector). If your data
suffers from an unknown spatially variable bias, this program will not
solve (all) your problems.
Directory structure
-------------------
The directory structure of the software is as follows:
README : this file
VERSION : text file that lists the software version
prog/ : source of the main program
sub/ : subroutines called by the main program
bin/ : executable of the main program
lib/ : object libraries
iraf/ : subroutines for reading and writing IRAF files
examp/ : examples of the use of the program
Installation
------------
To install the software edit the file setup.com. All site dependent
environment variables should be made to point at the correct directories
and libraries. Then issue the UNIX command
prompt> source install.com
This will construct an executable and put it in the directory bin/.
Once the programs have been installed and the executable created, the
environment variables in setup.com are not needed anymore. However, they
are required every time any of the programs is recompiled. Users who
will be modifying the programs regularly may wish to add the command
source /setup.com
to their .cshrc file (where should be the path to the directory
that contains the file setup.com).
Main Program
------------
There is only 1 main program, written in Fortran:
unpedestal.f : Reads a flat-field and a NICMOS image. Analyzes the two,
and returns both an estimate of the pedestal, and a new
image that has the best pedestal-value subtracted.
The executable of the program is called unpedestal.e, and is stored in
the directory bin/.
Examples
--------
The directory examp contains an example flatfield and NICMOS image (of
an uninteresting piece of sky obtained with the NIC2 Camera with
the F205W filter), as well as examples of the use of the program. To
run the example issue the UNIX command
prompt> source examples.com
Known Bugs:
-----------
Currently none.
*********************************
PART II: Discussion of the method
*********************************
Introduction
------------
After installation of the NICMOS instrument on HST it was found that
images suffer from an unknown, time-variable bias, which has come to
be known as the `pedestal'. A change in instrument operation on August
21, 1997 lessened the severity of the pedestal for all subsequent
data. Nonetheless, the pedestal is still present at some level in all
data taken with NICMOS. It should be noted that it is neither
particularly well established that the pedestal is indeed spatially
constant, nor that it is a pure bias (i.e., does not scale like the
flat-field). Nonetheless, the software discussed here operates on the
assumptions that the pedestal is indeed a constant bias.
Method
------
The program asks for the name of a calibrated (flat-fielded) data
image and the corresponding flat-field image. It then assumes that
there may be an unknown constant bias (the pedestal) in the
data, that was not removed in the calibration process. This means that
the calibrated image (C) equals the actual image (I), plus the
pedestal (P) times the flat field (F):
C = I + (P*F)
The program tries to determine the pedestal value. It exploits the
fact that the additional term (P*F) in the calibrated image introduces
a spread in the pixel values in the image. The program loops over a
range of trial values (T) for the pedestal, and for each it generates
a trial image I_T = C - (T*F). It applies an operation (O) to each
image to obtain an image J_T = O[I_T]. The operation is chosen to
remove unwanted features or spatial frequencies (see below). The
program finds the trial value T for which the spread of the pixel
values in J_T is minimal. The result is adopted as the best guess for
the pedestal.
The program provides various choices for the operation O, that defines how
J_T is obtained from I_T.
1: J_T = I_T
This implies that all Fourier frequencies in I_T
are used. This approach may work well if the image is free of
sources, for example, if it is a dark image.
2: J_T = a ring-median filtered version of I_T (i.e., every pixel
is replaced by the median value in a ring of pixels around it).
This implies that all pixel-to-pixel variations are removed,
and only the low Fourier frequencies in I_T are used.
This approach may work well for a field that is mostly empty, with a
few small sources (stars). This approach has been used by
Mark Dickinson in another package developed for pedestal removal.
3: J_T = I_T minus a median filtered version of I_T
This implies that low frequency structure is removed, and only
pixel-to-pixel variations remain. Thus only high frequency
information in I_T is used. This approach may work well
for any field, even if it contains stars or galaxies, or even it
contains a galaxy that fills the whole image.
To determine the spread S in the pixel values in a trial image J_T,
the program adopts the following approach. It makes an array of all
the pixel values. This array is sorted in ascending order. The modal
and median pixel value are determined, as well as the difference in
the pixel values corresponding to the percentiles 0.8413 and
0.1587. Half of this difference is used as an initial estimate
[gaudisp] for the Gaussian dispersion of the pixel values. An array is
then set up with Nbin bins between mode-(binmax*gaudisp) and
mode+(binmax*gaudisp). The pixel values are binned on this array to
obtain a pixel value histogram, and a Gaussian is fitted to this
histogram. The dispersion of the best-fitting Gaussian is used as the
spread S that we wish to minimize.
Using the above definition of the spread S, the program determines the
pedestal value that gives the minimum S using a straightforward
numerical minimization algorithm. (***) This result should be
interpreted with some care, since numerical errors (may) make the
function S(T) jagged near its minimum. The minimum found by the
routine may be a local minimum due to a `noise spike'. The program
therefore proceeds by mapping the function S(T). The value of S is
calculated on an array of (2*Nrefine)+1 values around the minimum,
using a stepsize pedstep. A smooth polynomial of order Morder is
fitted to these sample values, and the minimum of this polynomial is
determined. This minimum is adopted as the best guess for the
pedestal. The RMS deviation between the polynomial fit and the sample
values is calculated and used as an estimate of the numerical noise
S_err in S. The range within which the polynomial fit to S is smaller
than S_min + S_err is returned, and may serve as an estimate of the
numerical accuracy of the pedestal determination.
Bad pixels and outliers
-----------------------
It is important that stars, galaxies, cosmic rays or other sources are
removed from the analysis. The program uses two strategies to try to
to do this:
(a) It uses a median filter (as specified by the user) to remove
the most corrupted spatial frequencies.
(b) It minimizes a quantity S that is not very sensitive to
bad pixels and outliers, because:
(i) Gaussian fits are sensitive primarily to the
core of a distribution, and not to its wings;
(ii) pixels further than (binmax*gaudisp) from the mode are
not used.
This is a prime advantage over the alternative of using e.g.
the true second moment of the pixel value distribution as a
measure of the `spread'.
These properties make it less necessary for the user to generate a
pixel mask for each and every image that must be
analyzed. Nonetheless, if there are obvious sources that may corrupt
the analysis, it may not be a bad idea to mask these using a mask
image that can be given to the software on input. To identify whether
this is necessary, it is wise to first run the software without a
mask, and to then study the resulting *T.imh output images. These show
the result of applying the operation O to the input data and
flat-field images. Typically it is a good idea to mask the
coronographic hole, if present in your data (NIC2), as well as
severely vignetted regions of the image.
Removing differences between quadrants
--------------------------------------
After the best-fitting pedestal has been determined and subtracted, it
may be the case that there remain differences in the DC levels between
different quadrants. The program includes an algorithm that attempts
to remove these differences.
The user must specify a range ix0 to ix1 with ixo < ix1. These numbers
indicate the maximum row or column distance from a quadrant edge that
the algorithm will use in its analysis. For example, consider the
difference between the data to the left and to the right of the
vertical qudrant border in the middle of the image. If we set ix0=2
and ix1=10, then the program will compare the columns 119 to 127 (in
the quadrant left of the border) to the columns 130 to 138 (in the
quadrant right of the border). At each row of each of these two column
ranges the program will fit a polynomial to the data. The order of
this polynomial is Morcont, and can be specified by the user. The
polynomial fits on either side of the border are extrapolated onto the
border, and the difference between the extrapolated values from either
side is calculated. This is done for every pixel on the border. The
results for the 128 pixels along every border are then stored in an
array, and the mode of this array is taken as the estimate for the DC
offset between the quadrants on either side. This is done for each of
the 4 borders in the frame.
Notice that the choice Morcont=0 will simply compare the average pixel
value in the regions on either side of the quadrant edge. The program
provides two additional options that may be useful, e.g., in a crowded
stellar field, where it is important that outliers are not taken into
account. These options are to use either the median (Morcont=-1) or
the mode (Morcont=-2) of the regions on either side of the quadrant
edge.
The algortithm then searches for the DC offset p_i that must be added
to each quadrant i to best remove these differences. By definition the
algorithm searches for values such that the sum of the p_i is zero. So
there are only 3 free parameters, and 4 constraints. This problem is
solved in a least squares sence. Having solved this problem, the
algorithm provides two options:
a) add the calculated correction to each quadrant, and return.
b) add the calculated correction multiplied by the flatfield
to each quadrant, and return. This option is iterated 5 times,
to take into account the fact that the flatfield may not be
normalized exactly to unity in the region around the
quadrant boundaries.
***********************************************
PART III: Description of program in- and output
***********************************************
Input
-----
The program reads and writes in- and output images in the IRAF .imh
format. FITS files can be transformed to this format using the task
dataio.rfits in IRAF, or in IRAF 2.11 or later with the command
imcopy. The input images must have size 256 x 256 (as appropriate for
HST/NICMOS).
The input .imh files *must* contain numbers of datatype real (not
double precision). This will be OK if you simply copy a FITS file
produced by the pipeline to .imh format as described in the previous
paragraph. However, this may not be the case if your FITS files were
created by IDL or some other program. To see the datatype of your .imh
files, use the command `imheader long+' in IRAF. If you run the
software with .imh files that do not have the datatype real, you may
not get a meaningful error message; the program will just crash, or
produce nonsensical results. If you have FITS files with datatype
double precision that must be transformed to IRAF .imh files with
datatype real, then use the IRAF command dataio.rfits to do the
transformation.
The input data image *must* have been flat-fielded. The image may
either be an individual readout of an .ima image, or it may be the
.cal image in which all the individual NICMOS readouts have already
been combined with cosmic ray subtraction. Note that HST flat-fields
are always multiplicative.
The program asks for its input on the command line. No additional
files are required. The following information has to be supplied to
the program:
* How does the fit have to be done ? There are three options:
0: don't do a pedestal determination at all, i.e., fix the
pedestal to zero. This is useful for users that only wish to apply
the quadrant equalization algorithm.
1: determine the best-fitting pedestal for each quadrant separately.
The output image contains all four quadrants, but each quadrant
will have a different pedestal subtracted.
2: determine the best-fitting pedestal for all quadrants together.
The image is treated as though it doesn't have any individual
quadrants, and only one value of the pedestal is determined.
The output image contains all four quadrants, and each quadrant has
the same pedestal subtracted.
* In case of pedestal determination, the type of image for which the
program has to minimize the spread. 1, 2, or 3, as defined above.
* In case of pedestal determination, there is an option to
run the program in a quick-and-dirty manner.
In this case the program stops at the point marked (***) above.
This may not be appropriate for high accuracy work, but may be useful
for a quick look analysis.
* After the best pedestal has been determined (or fixed to zero)
and has been subtracted, there is
the option to apply an algorithm that attempts to make the transitions
between quadrants continuous
0 : Do not use this algorithm
1 : achieve continuity by adding a constant
to each quadrant
2 : achieve continuity by adding a constant
times the flatfield to each quadrant
* The user may supply a mask image. All pixels for which the mask image
has the value 1.0 will be used in the analysis. All pixels for which
the mask image has the value 0.0 will be ignored in the analysis.
The mask image should not contain any values other than 0.0 and 1.0.
* If the user wants to use non-default values for algorithm parameters,
he has to specify these. The following parameters may be specified:
{Parameters that are used by the entire program:}
[binf] : Fractional half-bin size for the mode finding algorithm.
To be able to calculate the mode of an array x(1), ..., x(N),
the program sorts the array, and then calculates the
density at each point using a running-average. The density
is defined as [binf * N] / dist, where dist is the binsize
that contains [binf * N] data points of the sorted array.
The default value for binf = 0.10. Larger values correspond
to bigger smoothing in calculating the density.
{Parameters that are used only by the pedestal determination part:}
[Nbin, binmax] : Used in the construction of a pixel histogram,
as defined above.
[ftol] : Requested fractional accuracy in minimization and root-finding
algorithms.
[Nrefine, pedstep]: values that define the array on which the function
S(T) is mapped, as defined above.
[Morder] : polynomial order used in fitting the sample values of S(T),
as defined above.
[MMin], [MMout]: The operation that is performed on the trial image,
as explained above, needs to calculate a median filtered image.
The median is calculated for a region between two squares.
Each square is defined as having x-values or y-values between
-M and M (inclusive). For the case of a normal median filter one
would set MMin = 0. For the case of a ring-filter one would say
MMin > 0. Note that the ring in this case is not a circular ring,
but a `square ring'.
{Parameters that are used only by the quadrant equalization part:}
[ix0], [ix1]: distance range of rows or columns from a quadrant edge
to be used in the algorithm that attempts to remove DC differences
between quadrants, as defined above.
[Morcont]: polynomial order to be used in the algorithm that
attempts to remove DC differences between quadrants, as defined
above. Morcont = -1 or -2 doesn't do a polynomial fit at all,
but uses the median or the mode of the pixel values, as
explained above.
Output
------
Images
------
The program generates a corrected output image. The output filename is
the input file name with either the letter A or the letter B
appended. The letter A indicates that the pedestal was determined for
each quadrant separately. The letter B indicates that only one
pedestal was determined for all quadrants simultaneously.
If the program is used for pedestal determination, then the it also
generates output images which have the additional letter T
appended. These contain the results of applying the operation O (as
defined above) to the input flat-field and the data image,
respectively. It is wise to visually examine these images for obvious
regions/sources that may have biased the results of the program. If
such regions and or sources are found, the program should be rerun
with these regions masked.
The program will not overwrite existing images. If it tries to create
an ouput image that already exists, it will generate the error message
`WRIFILFITS: error creating'
(but the program will continue running; it just won't overwrite the
output image). So when you wish to run the program twice in a a row,
make sure to first delete or move the output of the first run.
.out file
---------
The program generates an output file with the extension .out. This file
contains the same information that is written to the screen when the
program is running. It contains columns with the following
information:
[pedestal determination table]
output image : name of the output image
quad : quadrant number, defined such that images look as follows:
2 4
1 3
If only one pedestal is determined for all quadrants
simultaneously, then the number 0 is entered here.
mode : the mode of the pixel values in the input image.
median : the median of the pixel values in the input image.
pedestal : the best-fitting pedestal value, determined from a
polynomial fit to the function S(T) as described above.
In quick-and-dirty mode it is set equal to ped_bestD
as defined below.
new mode : the mode of the image after subtraction of the best-fitting
pedestal. Need not be exactly equal to mode - pedestal, because
the flat-field may not be normalized exactly to unity.
new median : the median of the image after subtraction of the best-fitting
pedestal. Need not be exactly equal to mode - pedestal, because
the flat-field may not be normalized exactly to unity.
ped_bestD : the best-fitting pedestal determined by a
numerical minimization on S(T), as described above. Should
be less accurate than `pedestal'.
spr_minD : the spread for the pedestal ped_bestD.
ped_lowP : lower bound on the pedestal, as obtained from the
polynomial fit, taking into account numerical noise as
described above. Set equal to the dummy value 0 in
quick-and-dirty mode.
ped_highP : upper bound on the pedestal, as obtained from the
polynomial fit, taking into account numerical noise as
described above. Set equal to the dummy value 0 in
quick-and-dirty mode.
spr_minP : the spread of the polynomial fit at its minimum. Set
equal to the dummy value 0 in quick-and-dirty mode.
Pfit_rms : the RMS residual of the polynomial fit. Set
equal to the dummy value 0 in quick-and-dirty mode.
Pfit_max : the maximum residual of the polynomial fit. Set
equal to the dummy value 0 in quick-and-dirty mode.
The region [ped_lowP,ped_highP] is always placed around pedestal, and
indicates the **numerical** error on the inferred value of pedestal.
However, it is not a true error estimate, since systematic errors may
play a role. The numerical error refers to the solution of the problem
that the program was requested to solve. This problem may not always
be relevant for the actual data; e.g., if the actual bias in the image
is not spatially constant, or if stars, galaxies or other sources have
not been properly removed by the median filtering operation.
[quadrant equalization table]
For each iteration there are two tables. The first has the following columns:
quad : quadrant number, as defined above.
correc : the calculated DC correction that has been added to the quadrant.
new mode : the mode of the pixels in the quadrant after the correction
has been applied.
new med : the median of the pixels in the quadrant after the correction
has been applied.
The second table has the following columns:
Quantity : Describes the quantity that is listed in the row.
In the first four rows this is the calculated difference
between two adjecent quadrants. In the last row
this is the RMS residual of the least-squares fit.
Goal : the calculated difference that the algorithm tries
to remove.
Fit : the actual fitted difference that the program has been
able to remove by using the corrections listed above.
Goal_med : Similar as goal, but now taken to be the median of the
array of extrapolated differences along the border, rather
than the mode (as in the definition of Goal). These values
are not used by the program, and are returned only to signal
possible problems (if Goal and Goal_med are very different,
then probably something is wrong).
.log file
---------
If the program is used for pedestal determination, then it also writes
a .log file (but not in quick-and-dirty mode). This file contains
calculations of the spread as function of the pedestal value. These
can be plotted at a later time with any regular plotting package, to
visualize what the program has done. The following columns are listed
in this file:
quad : quadrant number. If only one pedestal is determined
for all quadrants simultaneously, then the number 0 is entered
here.
ped_trial : the trial value for the pedestal.
spread : the calculated spread, as defined above.
spread_fit : the polynomial fit to the spread.
residual : spread - spread_fit.
Error messages -- Array sizes
-----------------------------
The parameters listed in paramval.h in the directory /prog determine
the sizes of the arrays that are reserved by the programs. When the
input data is too large for these arrays an error message is
generated, and the program aborts. This should not happen, because
HST/NICMOS images always have the same size (256*256).
************************************
PART IV: Experience with the program
************************************
Here are some preliminary results obtained by Roeland van der Marel &
David Zurek. We focused on Method 2 (using the operation number 2
defined above, i.e., J_T = a ring-median filtered version of I_T) and
Method 3 (using the operation number 3 defined above, i.e., J_T = I_T
minus a median filtered version of I_T). Method 1 is less likely to be
useful, unless in special cases.
a. Methods 2 and 3 give consistent results when applied to an image
with no sources at all.
b. Method 3 works equally well when there no sources at all, when there are
stars in the image, or when a galaxy fills the whole image.
c. Method 2 does not work if a galaxy fills the image, or if any other
low-spatial frequency variations are present in the image. It is unknown
whether ther are any situations/images for which Method 2 is more
accurate than Method 3.
d. With the above methods, the pedestal can be determined from an
individual image to an accuracy of approximately 0.01-0.02 DN/sec RMS.
e. The pedestal does not vary among expsoures in a given visit by more
than 0.01 DN/sec RMS (for the same camera, filter, read-out sequence,
etc.)
f. The absolute pedestal level changes with time of year. Part of this may
be due to time variability of the shading.
g. The relative pedestal level from quadrant to quadrant doesn't change
(much) with time of the year.
*******************
PART V: Suggestions
*******************
The program has many switches that the user can play with to optimize
the results for the particular application at hand. However, not
everyone will want to invest time in this, and for some applications
only modest accuracy will suffice. For such situations I recommend the
following settings (see the discussion under `input' above), which
appear to give good performance in many cases that I have tried :
* how to do the fit: 1
* type of image for which to minimize the spread: 3
* quick-and-dirty mode: 0
* quadrant equalization: 1
* use default parameters: 1
These options are indicated by the program when it runs, as `NOVICE
USER HINTS'.
If you are not satisfied with the results of the program when run with
these options, you can proceed by making changes in the default
settings. The setting of the parameters in the quadrant equalization
algorithm is probably most critical. The optimum settings will depend
quite a bit on the structure of the sources in your image. If the
default settings do no make the quadrants as continuous as you would
like, then familiarize yourself with the quadrant equalization
algorithm (see above), and try non-default values for ix0, ix1 and
Morcont. If the quadrant equalization doesn't work properly in a field
with no extended sources, but with many point-like sources, it may be
useful to try the non-default value Morcont = -1 or -2. Note that
Morcont = -2 will compute and use the mode of ix1-ix0+1 pixels, which
is not very well defined if ix1-ix0+1 is small (less than 10-20 or
so).
****************
PART VI: History
****************
Here is a brief history of the versions of this software:
1 (Jan-Mar 1998) : Initial release version
2 (April 1998) : Added quadrant equalization algorithm
3 (May 17, 1998) : Added use of Morcont = -1, -2 in the quadrant equalization
algorithm (see README for details).
4 (May 20, 1998) : Added option to skip pedestal determination, and
apply only the quadrant equalization algorithm.
5 (July 27, 1998): Minor change to avoid generation of an error message
that occurred in certain pathological circumstances.
Roeland van der Marel.
Baltimore, July 27, 1998
Return to my home page.
Last modified December 8, 1998.
Roeland van der Marel,
marel@stsci.edu.
Copyright © 1998 The Association of Universities for
Research in Astronomy, Inc. All Rights Reserved.