Subject: The meaning of variance maps



Dear Roland:

Thanks for your mail on the mosaic imaging issue. It has taken me
a while to answer because things I believed to understand were not
immediately so obvious as I thought.

We discussed my weighting scheme here last week and Niels-Joergen 
suggested a simplification which avoids the issue of the order of 
the addition process:

Niels-Joergen argues that in forming the weight of each image
one does not need to consider the exposure already accumulated
during the previous steps. His point is that for each pixel in
the final image there will be a certain total exposure and all
contributions to this pixel should be weighted relative to this
accumulated exposure. This normalization can be done after all
contributions have been included because it is part of the
existing procedure that we maintain a map of the accumulated 
exposure for each pixel. When the weights only depend on the
data of each image by itself it is obviously unimportant in
which order the addition is done (this is, I believe, what you
call "symmetry" in your mail).

*** Midisky variance maps:

Consider the uncertainties in the sky image from a single observation:

1. For a single observation we know how much effective time, T,
(including dead time and grey filter effects) we have spent obtaining 
N_i photons backprojected to the i-th sky pixel. In the simple 
backprojection scheme all the photons arriving at a given sky pixel 
are independent and therefore should obey Poisson statistics. However, 
(almost) all sky pixels have some photons in common and are therefore 
not independent measurements.

2. We also know how much detector area, A_i, actually was illuminating 
this pixel (from the illumination map). The product of T and A_i is
the exposure E_i associated with the N_i observed photons in pixel i.

3. From the total number of counts, B_i, in a set of k neighboring sky
pixels we may obtain an estimate of the background in the i-th pixel.

Thus the excess counts in the i-th sky pixel is estimated as:

   S_i = N_i - B_i/k                                (1)

with a standard deviation:

   sigma_i = sqrt(N_i + (B_i/k)^2)                  (2)
   
(Poisson statistics and large numbers assumed).


(As I see it it does not really matter which method is used to
subtract the background - balanced reconstruction or subtraction
of the locally estimated background. Both methods have their
pitfalls if applied without care. In midisky I use the latter,
I select 20 pixels in a fixed pattern around the i-th pixel 
and make an average over 18 of these, discarding the two pixels
having the highest counts - avoiding to create negative artefacts
around strong sources).

If we want to normalize our image in counts/s/cm2 we must apply
the exposure factor to both S_i and sigma_i:

   I_i:    (S_i +/- sigma_i) / (A_i * T)            (3)

*** Mosaics

The weighting algorithm I proposed should form a mosaic image with 
the maximum signal to noise in each pixel. This should be optimum
for source detection purposes. But this implies that the weighting 
of each image depends on the variance in the image, and this again 
depends strongly of the number and strength of sources within the 
FOV. So the weight factors for a given source depends on the activity 
level of the surrounding sources. Therefore it is not possible to 
derive accurate fluxes from these images (neither absolute nor 
relative).  

Nevertheless, the summation formulaes I have given for signal, 
exposure and variance should maintain the relation between variance 
and signal so the formulaes (1) to (3) should be applicable to the 
final image. The number of counts in the variance maps of a mosaic 
image is the weighted sum of the counts in the individual images.
However, and this is where I had difficulties: With Niels Joergens new
proposal for the practical weighted addition, I have difficulties
arguing for the precise meaning of of the summed variance maps.
Below I offer my present understanding.

Niels Joergen proposes to form the sums for each pixel in the
mosaic (for pixel index i and input image index j): 

S_mosaic_i = SUM_j{w_ji * S_ji}                     (4) SIGNAL
E_mosaic_i = SUM_j{w_ji * E_ji}                     (5) EXPOSURE 
V_mosaic_i = SUM_j{w_ji * w_ji * V_ji}              (6) VARIANCE

with weights: w_ji = E_ji / V_ji.                   (7) WEIGHT

In order to be able to complete the normalization after the 
mosaic image is formed is is now also necessary to maintain 
a fourth and fifth image containin the sum of the weights
and the sum of the weights squared:

w_mosaic_i  = SUM_j{w_ji}                           (8)
w2_mosaic_i = SUM_j{w_ji * w_ji}                    (9)

Although I am not able to describe the meaning of the mosaic 
variance map in the general case, at least it is simple to see
that the above formulae gives a sensible result if you add M
identical maps, then:

S_mosaic_i  = M * w_ji * S_ji                      (10  
E_mosaic_i  = M * w_ji * E_ji                      (11)
V_mosaic_i  = M * w_ji * w_ji * V_ji               (12)
w_mosaic_i  = M * w_ji                             (13)
w2_mosaic_i = M * w_ji * w_ji                      (14)

and the end result (after dividing each sum with the mean value 
of the corresponding weights) becomes:

S_mosaic_i = M * S_ji                              (15)
E_mosaic_i = M * E_ji                              (16)
V_mosaic_i = M * V_ji                              (17)

So the signal/noise improves as the square root of M.

The division with the average weight (and not with the sum of the
weights) is required if you want the values in the mosaic maps to
to increase when you add more images. If we divide with the sum
of the weights we will arrive at the average single sky image
contribution - which is not what we want. If we now assumes that
this procedure also applies in the general case, then we have a 
full prescription for generating mosaics, and I will argue that
the final variance maps represent a good measure of the effective
number of photons forming the signal S_mosaic_i. And I will also 
argue that S_mosaic_i expresses the integrated effective excess 
number of photons in pixel i, with an estimated statistical
error (sigma) of sqrt(V_mosaic_i).

If you plot S_mosaic by itself you will get an image where the 
signals are large in the center and small at the edge. This is 
actually a very nice mage to look at, but you cannot at all 
compare the strengths of the different excesses because you have 
not accounted for the differences in exposure. But I believe this 
image is actually the best in which to spot new sources.

If you plot S_mosaic/E_mosaic you compensate for the differences
in exposure for the various parts of the image, but now you
blow up the the small signals (with their correspondingly large
relative fluctuations) near the edge, and the image becomes
much less nice to inspect. My proposal for practical work is
 to work with S_mosaic when you search for sources, and with 
S_mosaic/E_mosaic for presentation work but then limit the
displayed area to the central parts of the mosaic where the 
exposure is good.

with best regards

niels