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 axes12.
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-variances13. 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:
where is the expected intensity in image , is the analyser transmission factor, is the analyser efficiency factor, and 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 . The Stokes parameters , and are chosen to minimize , half the weighted sum of the squared residuals between the expected intensity and observed intensity , where the sum is taken over all images (i.e. all values of ). The weights used are the reciprocals of the variances associated with the observed intensity values:
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 multiple14 of the standard deviation associated with the observed intensity value, , 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 for (i.e. give all input values equal weight). In this case no variances for , and 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:
Each step in this process is described further below:
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.
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.
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, is larger than 3, then the residuals produced in 7 will already represent in some way the mean of noise samples, and so the need to smooth the squared residuals is less.
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.
12The only exception is that when data is “smoothed” it is only smoothed in the two spatial axes - no smoothing occurs between spectral channels.
13Neither do they use the improved de-biassing technique described by Sparks & Axon.
14Specified by POLCAL parameter NSIGMA.