E.1 The Noise Level in Individual Input Images

E.2 Using POLSIM to Investigate Noise Characteristics

E.2 Using POLSIM to Investigate Noise Characteristics

This section describes the algorithm used by POLCAL to created Stokes vectors from single-beam
data. It assumes that the data is single-waveband. Spectropolarimetry data is processed in the same
way, except that an intensity ‘image” is actually an intensity “cube” containing two spatial axes and
one spectral axis. In this context, no distinction is usually made between spatial and spectral
axes^{12}.

Given a set of observed intensity values and corresponding variance values, a set of Stokes vectors
with associated variances can be produced following the technique described by Sparks & Axon
(*P.A.S.P.* in prep.). The full error analysis given by Sparks & Axon also allows the estimation of the
co-variance between Q and U which occurs when non-orthogonal analyser positions are
used. POLCAL produces these co-variances and stores them as a separate 2D NDF within
the POLPACK extension of the output cube. However, the calculation of the degree and
orientation of the polarization performed by POLVEC and POLBIN do not as yet use these
co-variances^{13}.
It is hoped that this will be rectified in a future release of POLPACK.

The technique described by Sparks & Axon is basically a case of fitting the intensity values at each pixel with the following function:

$$\begin{array}{rcll}{I}_{k}^{\prime}& =& \frac{t}{2}\left(I+\u03f5.\left(Q.cos2{\varphi}_{k}+U.sin2{\varphi}_{k}\right)\right)& \text{(2)}\text{}\text{}\end{array}$$

where ${I}_{k}^{\prime}$ is the expected
intensity in image $k$,
$t$ is the analyser transmission
factor, $\u03f5$ is the analyser
efficiency factor, and ${\varphi}_{k}$
is the anti-clockwise angle from the reference direction to the analyser
position (or effective analyser position if using a half-wave plate) for image
$k$. The Stokes
parameters $I$,
$Q$ and
$U$ are chosen to
minimize ${\chi}^{2}$,
half the weighted sum of the squared residuals between the expected intensity
${I}_{k}^{\prime}$ and observed intensity
${I}_{k}$, where the sum is taken over
all images (*i.e.* all values of $k$).
The weights used are the reciprocals of the variances associated with the observed intensity
values:

$$\begin{array}{rcll}{\chi}^{2}& =& \frac{1}{2}\sum _{k}\frac{{\left({I}_{k}-{I}_{k}^{\prime}\right)}^{2}}{{\sigma}_{k}^{2}}& \text{(3)}\text{}\text{}\end{array}$$

A minimum of three images are required, taken at different analyser positions, to perform this fit. The
basic algorithm gives little protection against badly aberrant input values (*e.g.* cosmic rays, etc.)
which can corrupt the fit. However, if more than three images are obtained with different
half-wave plate positions, then an iterative process can be used to identify and reject such
aberrant data values. This process works by first calculating the Stokes vectors using all
input data values. The expected intensity value at each input pixel is then found from
these Stokes vectors using equation 2. The residual between this expected intensity value
and the observed intensity value is then found. If this residual is larger than a specified
multiple^{14}
of the standard deviation associated with the observed intensity value,
${\sigma}_{k}$, then
the observed intensity value is rejected. The Stokes vectors are calculated again, excluding all
the rejected intensity values. This process is repeated a number of times, as specified by
POLCAL parameters MAXIT and TOLR. The number of pixels rejected in each image at
each iteration can be displayed by setting POLCAL parameter ILEVEL to 3. There is an
option to reject output pixels (*i.e.* set them bad) if more than a given fraction of the input
intensity values are rejected as aberrant by the above iterative procedure (see parameter
MINFRAC).

If no usable variances are available for the observed data values, the simplest option is to assume a constant
value of $1.0$
for ${\sigma}_{k}$
(*i.e.* give all input values equal weight). In this case no variances for
$I$,
$Q$ and
$U$ can be
produced, and the iterative rejection of input data described above cannot be performed. This option
can be used by setting POLCAL parameter WEIGHTS to 4.

In the absence of input variances, another option is to estimate these variances by looking at the deviations of the observed intensity values from the best fitting sine curve given by equation 2. This allows variances for the output Stokes vectors to be created, and also allows the iterative rejection of bad data. This option can be used by setting POLCAL parameter WEIGHTS to 3. The technique is described in the following pseudo-code, and is illustrated diagrammatically in Figure 11:

1 Fill the variance image with the value 1.0

2 Calculate Stokes vectors

3 Set ITER = 0

4 While ( not converged and ITER < MAXIT ) do...

5 Reduce the noise in the Stokes vectors

6 Calculate expected intensity values implied by the smoothed

I, Q and U values

7 Find the squared residuals between these intensity estimates

and the observed intensity values

8 Find the mean squared residual at each pixel position

9 Smooth the image holding these mean squared residual values

10 Copy this smoothed image to the current variance image

11 Copy the input intensity values, excluding values for which

the squared residual found above is too large

12 Calculate I, Q and U again using only the copied input values,

and the current variance image

13 Set ITER = ITER + 1

14 End do

15 Reject any Stokes vectors corresponding to pixel with few good

input intensity values.

2 Calculate Stokes vectors

3 Set ITER = 0

4 While ( not converged and ITER < MAXIT ) do...

5 Reduce the noise in the Stokes vectors

6 Calculate expected intensity values implied by the smoothed

I, Q and U values

7 Find the squared residuals between these intensity estimates

and the observed intensity values

8 Find the mean squared residual at each pixel position

9 Smooth the image holding these mean squared residual values

10 Copy this smoothed image to the current variance image

11 Copy the input intensity values, excluding values for which

the squared residual found above is too large

12 Calculate I, Q and U again using only the copied input values,

and the current variance image

13 Set ITER = ITER + 1

14 End do

15 Reject any Stokes vectors corresponding to pixel with few good

input intensity values.

Each step in this process is described further below:

- (1)
- The
*variance image*holds the variance to associate with each input intensity value. Initially, all variances are set to $1.0$ so that all input data values have equal weight. Note, there is only one variance image. A given input pixel is given the same variance in equation 3, no matter which input image it comes from. If this were not the case (*i.e.*if each input image had its own variance image) the better images would receive higher weight, resulting in them being fitted more closely. On the next iteration, they would then be found to have smaller residuals, and so would receive even higher weight, resulting in them being fitted even more closely. This circular process would repeat until the fit was totally dominated by the best input image. - (2)
- The Stokes vectors are calculated using the constant weights set up in step 1. All input data is used.
- (3)
- The variable ITER counts the number of iterations which have been performed.
- (4)
- Loop round performing iterations until the process converges, or the maximum number
of iterations specified by parameter MAXIT is reached. If MAXIT is zero, the Stokes
vectors produced in step 2 are returned.
The convergence criterion is specified by parameter TOLR. Convergence is assumed to have been reached when the number of pixels rejected from each input image changes by less than TOLR pixels between iterations. If the number of pixels rejected from any image changes by more than the value of TOLR, then another iteration is performed (subject to the MAXIT limit). The default value for TOLR is zero, resulting in iterations being performed until the number of pixels rejected from each image no longer changes.

- (5)
- The three images holding $I$,
$Q$
and $U$
are smoothed to reduce the noise. The extent of the smoothing (in pixels) is controlled
by POLCAL parameter SMBOX. Setting SMBOX to 1 or 0 results in no smoothing being
performed.
As mentioned above, the likely error in a set of observed intensity values can be estimated by comparing them with the best fitting sine curve defined by the current Stokes vectors. This estimation can be performed independently for each pixel. However, unless your images are badly under-sampled, it is also reasonable in the absence of noise to expect each image to be spatially smooth over some small scale. So we could also estimate the noise by spatially smoothing each image and looking at the residuals between the smoothed and original data. The two methods can be combined by smoothing the Stokes vectors before finding the estimated intensity values.

In some cases this smoothing is essential. For instance, consider the case where only three analyser positions are available. In the absence of any smoothing, the best fitting sine curve would always be an exact fit to the three observed data values, no matter what the noise may be. Therefore the input variances could not be estimated. Introducing some spatial smoothing changes the expected intensity values so that they are no longer identical to the observed values, thus producing non-zero residuals and allowing the input variances to be estimated. In general, the fewer the number of analyser positions, the more important is the spatial smoothing.

However, spatial smoothing can introduce problems. The loss in resolution can result in real structure being interpreted as noise (

*i.e.*having large residuals). This in turn, can result in the input variances being over-estimated, and pixels being rejected from the calculation of the Stokes vectors. This is particularly apparent close to small bright features such as stars, but generally affects any region with significant spatial gradient. For this reason, you may want to consider turning off the smoothing if you have sufficient analyser positions to manage without it. Having said that, the smoothing is implemented in a way which minimises these effects by minimising the loss of resolution. A common method for smoothing an image is to convolve the image with some Point Spread Function. For instance a block filter replaces each pixel with the mean of the values within a box centred on the pixel. The smoothing performed by POLCAL is not done in this way. Instead, a quadratic surface is fitted to the data in the box, and the value of the fitted surface at the central pixel is taken as the smoothed pixel value. This method is illustrated in Figure 12.

- (6)
- Equation 2 is used to calculate the intensity values which would be expected on the basis of the smoothed Stokes vectors.
- (7)
- For each input image, the residuals between these expected intensity values and the corresponding observed intensity values are found, and squared. If the input images have differing zero points (caused by instrumental effects or imperfect sky subtraction, for instance), then the residuals for an individual image may well not have a mean of zero. If left uncorrected, this would result in the mean squared residual at each pixel being an over-estimate of the actual input variance. To reduce this effect, there is an option (see parameter DEZERO) to apply a zero point correction to each input image. This takes the form of a constant value to be added to the image. This constant value is chosen so that the mean of the residuals calculated in this step is zero, and is displayed if parameter ILEVEL is set to 3 or more. If the DEZERO option is selected, the zero point corrections are also used when calculating the Stokes vectors themselves. The DEZERO option slows down the convergence of the iterative procedure, resulting in a larger number of iterations being required, and thus a larger run time.
- (8)
- The images containing the squared residuals are stacked together into a single image holding the mean squared residual at each input pixel.
- (9)
- The image produced by the previous step is smoothed by replacing each pixel with the mean of
the values in a box centred on the pixel. The size of the box used is the same as in step 5. So
setting parameter SMBOX to 1 or 0 will result in no smoothing being performed at this
point.
Consider again the case of three analyser positions. The residuals found in step 7 will represent samples from the noise distribution. The mean of the squared residuals in a small box thus gives an estimate of the variance of the noise distribution. If the number of analyser positions, $N$ is larger than 3, then the residuals produced in 7 will already represent in some way the mean of $N$ noise samples, and so the need to smooth the squared residuals is less.

- (10)
- These mean squared residuals are used as the variance estimates for the next iteration. Note again, the same variance value is used for a given pixel, no matter which input image the pixel comes from.
- (11)
- The residual between the observed and expected intensity values at each pixel is compared with the square root of the current variance estimate for the pixel. Any pixel for which the residual is greater than NSIGMA times the standard deviation is flagged as bad, and is not used in the next iteration.
- (12)
- A new set of Stokes vectors are produced using the observed intensity values, and the current estimate of the variances. Any values flagged as bad in the previous step are not included in the calculation.
- (13)
- Indicate that another iteration has been performed.
- (14)
- Go round for the next iteration.
- (15)
- Once convergence (or the maximum number of iterations) has been reached, a check is made on the number of good intensity values which contributed to each output Stokes vector. If too many of the input intensity values were rejected as aberrant in step 11, then the Stokes vector is set bad. This criterion is set by parameter MINFRAC, which gives the minimum fraction of the input intensity values which must pass the check in step 11 in order to create a good Stokes vector. The default value of zero result in no Stokes vectors being rejected.

If the ILEVEL parameter is set to 3, the number of pixels rejected from each input image at each iteration is displayed. In addition, the RMS value of the residuals in each input image is displayed. This gives an estimate of the mean noise level in each input image. Note, though, that these RMS values are not used in the above algorithm, which gives equal weight to all input images.

There is an option for POLCAL to store these mean variance values in the VARIANCE components of the input images (see parameter SETVAR). If this option is selected, the VARIANCE array in each input image is filled with a constant value equal to the mean variance in the image estimated on the final iteration.

So you may want to run POLCAL twice; the first time just to estimate the mean noise level in each input image, and the second time to use these mean noise levels to calculate the Stokes vectors. The first time you set parameter WEIGHTS to 3 and SETVAR to TRUE, causing the input variances to be estimated and stored in the input images. You then re-run POLCAL setting parameter WEIGHTS to 1 in order to use these variances. This would then produce Stokes vectors in which the better input images have higher weight.

The POLSIM application is a useful tool for investigating the noise within a set of observed intensity images. It takes in a cube of Stokes vectors as produced by POLCAL, together with a set of template intensity images, and creates a set of intensity images analogous to the templates, but filled with the expected intensity values given by equation 2. The analyser properties (angle, efficiency and transmission) are defined by the templates, as are the required pixel positions.

So if you use POLCAL to generate a cube of Stokes vectors from a set of observed intensity images, you could then use POLSIM to convert the Stokes vectors back into intensity estimates, using the original observed intensity images as templates. In a noise-free world, the resulting estimated intensity images should be identical to the original observed data. The presence of noise in the input data will result in this not being the case, and the residuals between the two sets of images will give an indication of the noise in the input images.

POLSIM can also be used to verify that the Stokes vectors produced by POLCAL are consistent with the observed data.

^{12}The only exception is that when data is “smoothed” it is only smoothed in the two spatial axes - no smoothing occurs
between spectral channels.

^{13}Neither do they use the improved de-biassing technique described by Sparks & Axon.

^{14}Specified by POLCAL parameter NSIGMA.