General Notes and Advice for Fitting Images

WARNING: this page is currently an incomplete work in progress!

Pixel values and uncertainties

It is important to have some idea of the underlying statistical model for your data.

The basic assumption is that the pixel values in a data image are the outcome of some statistical sampling process (usually Poisson or Gaussian) operating on an underlying model of the intensity distribution.

A very simple example of an underlying model might be a constant sky background plus an elliptical 2D Sersic function representing a galaxy. Sampling of this by a 2D detector results in pixel values which are a combination of the image model and a sampling model for the distribution of individual intensity value recorded for each pixel.

If your data processing pipeline has produced some kind of error image (sigma or variance values) which you trust, then you can go ahead and tell imfit about the error image (via the --noise command-line option, plus the --errors-are-variances flag if the pixel values in the error image are variances instead of sigmas).

For the common case where you do not have an error image, you can have imfit estimate the per-pixel uncertainties for you, either from the data image or from the model. The question is then whether to use proper Poisson statistics or the usual Gaussian approximation of Poisson statistics (i.e., sigma ~ sqrt(intensity).

If the original count levels (including background) are high (say, greater than 100 photo-electrons per pixel), then you can probably assume the Gaussian approximation of Poisson statistics is OK, and use some variant of chi^2 as the fit statistic.

IMPORTANT NOTE: Converting your pixel values to photo-electrons (or particles, or…)

In order for imfit to correctly estimate the per-pixel uncertainties, it is very important that the pixel values be convertible to units of the original detected quantities – i.e., integrated photo-electrons for a typical detector, or particles/pixel for an N-body simulation.

Since the actual recording of an astronomical image almost always involves conversion of detected photo-electrons to counts (also called “data numbers” or ADUs) via an A/D converter, the output image is already no longer in the proper units.

So you must tell imfit how to convert the pixel values back to the original values, at least in a general, image-wide sense. (It’s probably overkill to worry about the effects of flat-fielding.)

The simple way to do this is to give imfit an effective gain factor, via the GAIN parameter in the config file, or the --gain command-line option. The effective gain is whatever number will, when multiplied by the image’s pixel values, convert them back to, e.g., photo-electrons. For a processed (but not photometrically calibrated) image, where the pixel values are in ADUs, this is simply the gain of the original A/D conversion, in units of electrons/ADU. Imfit will then multiply the pixel values by the GAIN parameter to turn them back into photo-electrons. (Note that some data-taking setups will record an inverse “gain” in units of ADU/electrons; you will need to invert this before passing it to imfit!)

If your image has photometrically calibrated flux units (Jy, nanomaggies, whatever), then the GAIN parameter must be a number that will convert these back to photo-electrons. If an image happens to be in units of magnitudes per square arc second, then it must be converted to linear intensity values first.

Special Notes for SDSS Images

DR7 (and earlier) Images

DR7 images are fairly straightforward; the only real potential for confusion is with the background level.

DR7 images do not have background subtraction applied (a crude estimate of the background level is included in the header of each image). However, each image has an artifical “soft bias” level of 1000 added to each pixel. Thus, if an image has a mean background level of 1210 counts/pixel, the actual observed sky background is only 210 counts/pixel. The additional constant value of 1000 should be removed and not included in any sky-background levels input to imfit.

Thus, if you determine that a particular image has a mean background level of, say, 1210.7 counts/pixel, you can either:

  • Subtract 1210.7 from the image, and use 201.7 (not 1210.7) as the ORIGINAL_SKY (or --sky commandline option) parameter for imfit;

  • OR: Subtract 1000 from the image, then include a background component (e.g., FlatSky) in the model with initial value = 210.7 (possibly fixed to that value, if you don’t want imfit to vary the background).

The A/D gain values are included in the tsField FITS tables that go with each field; typical values can also be found in the table at the bottom of this page

DR 8+ Images

Images that are part of DR8 and later releases are more problematic, because they have been preprocessed in slightly complicated and potentially confusing ways.

First, the pixel values have been photometrically calibrated and transformed into flux values – specifically, “nanomaggies” (units of 3.631e-6 Jy). To fit the images properly, you will need to convert these back to counts, or else combine the flux conversion with the A/D gain to make an “effective gain” parameter as input to imfit.

The photometric calibration is recorded in the image header using the NMGY header keyword. You can convert the image back to ADUs via

image_ss_counts = image_ss_nmgy / NMGY

Second, a 2D model of the background is computed and subtracted from the image as part of the pipeline, so you will definitely need to include some approximation of the background value as an ORIGINAL_SKY parameter in the configuration file (or via the --sky command-line option). The 2D background model is included in the image files in a somewhat obscure fashion, as a series of floating-point values for interpolation, stored as a table in extension 2 of the An additional problem is that these sky-model values are in counts, not in the nanomaggy units of the processed data image.

An estimate of the original sky background can be obtained by averaging the interpolation values in the second extension. For example, using numpy and the module in Python:

>>> import numpy as np
>>> from import fits
>>> hdu_list ='<path-to-SDSS-image-file>')
>>> sky_bintable = hdu_list[2]
>>> np.mean(['ALLSKY'])

The simplest approach is probably to convert the image from nanomaggies to counts, then use the standard A/D gain values and the mean sky value from the second FITS extension.

Alternately, you could combine the NMGY value with the gain to get an effective-gain value that converts the nanomaggie values directly to photo-electrons (“gain” = A/D gain / NMGY) – but then you will have to convert the sky value from counts to nanomaggies.