The Code

The algorithm, idac, has been written in c and is command line driven. An example of the command line is:

(Note that ALL files are in FITS format.)
 
 
conv
 This is the input convolution image data cube, i.e. N x N x M where N is the side of the array which should be in powers of 2 and M is the number of frames in the data cube. This is the observation after preprocessing, i.e. sky substraction, dark substraction and flat-fielding.
object
 This is the initial input object estimate. For the first start-up of the algorithm a mean of the observations usually represents a good start. This can be either a straight mean across the data cube, a centroided mean, or a peak-tracked mean. The dimensions are N x N.
 
psf
 This is the initial PSF estimates. If the initial estimate for each frame of the data cube is the same then this can be inputted as a single frame of size N x N. However, if the initial estimate differs for each frame of the data cube then this should also be a data cube of size N x N x M.
 
Support Constraints
 The algorithm will also expect to see three additional files at this stage. These are the convolution image, object and PSF supports and take the name of their respective input files appended by "_sup", i.e. conv_sup, object_sup, and psf_sup. If the supports are constant over the data set, then these can be set as singles frames of size N x N and are binary files with values of unity indicating that the pixel will be used in the fidelity constraint and zero that it will not. If the support changes from frame-to-frame, then conv_sup can be set to a data cube of size N x N x M. conv_sup can be used to describe a bad-pixel mask for image domain reductions.
 
weight
 This is an image of size N x N [x M] which gives the weighting of each pixel in the observation. An SNR weighting can be used here which directly modifies the convolution constraint. This is applied in image domain application of the code such that the pixels represent the weight at each pixel.
iter
 This is the maximum number of iterations for the algorithm to run.
 
mod
 This number determines the frequency with which the output files are updated. This is a useful parameter. For example, if it didn't exist then there would be no output if there were a system crash. Typically the user should set it to be about 10. However, this can vary and experience using the algorithm on different platforms can better determine its value, based on the algorithm's speed.
 
control
 This is probably the must important parameter and care needs to be taken in setting it correctly. This is a decimal number whose value is determined by the different options available for the algorithm. Each option sets a different bit. The following table shows which bit controls which option. When using a combination of options, simply add the decimal values together. A more thorough description of each option is given below.
 

BIT Value	Decimal Value	Constraint
bit 0	1	PSF Band Limit
bit 1	2	Mean PSF
bit 2	4	Not implemented
bit 3	8	Fixed object or PSF (w/bit 0)
bit 4	16	Sky Fitting  (Sky is ALWAYS subtracted)

PSF Band Limit - bit 0

 This flag determines if the band-limit constraint is used with the PSF. This constraint penalizes the solution if there are significant values at spatial frequencies of the PSF greater than the telescope diffraction-limit. It is recommended that this flag be set at all times as it prevents the PSF from converging to a delta function. An elliptical PSF is allowed but currently only if it is aligned along the x and y-axes. The cut-off frequency values have to be set in the FITS header of conv. The keywords are R0X and R0Y and typically are set the the same. The following is an example of a FITS header having these values set where these define the radius of the cut-off-frequency in pixels.  PA allows the position angle of an elliptical support region to be specified.

R0X = 64
R0Y = 64
PA   = 0

These values can be determined from the mean Fourier modulus (the square-root of the power spectrum) of the convolution images or measured PSF's. Note that the power spectrum can be easily computed using the powerspec tool in the STSDAS fourier package.

Mean PSF - bit 1

 This option penalizes the solution if the mean PSF estimate departs from a mean value. This mean value could be an estimate of the PSF obtained from a reference source or even a star in the field. It could even be a model of the PSF. When this flag is set, the algorithm expects to see a file with the same name as the initial PSF estimate but with the extension "_saa", i.e. psf_saa. This option is useful to use when the user believes that they have a good estimate of the PSF but not good enough for a standard known PSF deconvolution. The algorithm can be allowed to reach convergence with this flag set and then can be restarted with it switched off to permit further departures of the individual PSF estimates. psf_saa is of size N x N

The contribution of the Mean PSF constraint to the total error metric can be weighted between 0 and 1 depending upon the user's confidence in the mean PSF. A value of 0 indicates no confidence and mean PSF error metric will have no effect on the total error metric. A value of unity indicates full confidence. This weighting term appears as the keyword ALPHA in the FITS header of conv, the default value is ALPHA = 1.

Fixed Object or PSF - bit 3

 This option permits the algorithm to either fix the object or the PSF to perform a standard deconvolution. As to which is fixed depends on whether the PSF Band limit is also set. To fix the PSF, both the PSF Band limit and Fixed options must be set, i.e. control = 9. To fix the object, the PSF Band limit option must be off so that control = 8.

Sky Option - bit 4

 The algorithm can also estimate a mean sky background for each of the input convolution frames. This is especially useful when dealing with very noisy data as the inherent non-negativity constraint of the algorithm produces a half-wave rectification of the noise. When this flag is set, the algorithm uses the FITS image of size NxN called sky having a value equal to estimated the value of the mean sky for the data. An observed sky (flat-field) could be appropriate for this. 

Note: If bit 4 is not set, sky is simply subtracted (no fitting occurs). A sky file (sky) with all pixels set to zero must ALWAYS be provided when no sky fitting is required, as sky is always subtracted.
 

newstart
 This parameter is ignored in V2.7. However a numeric value still needs to be used, typically use "0".
begin
 The is an integer number representing the first frame of a range of frames within the data cube the user wishes to reduce. Note that as this is a c routine then the first frame is frame number 0.
end
 The is an integer number representing the last frame of a range of frames within the data cube the user wishes to reduce. Note that as this is a c routine then the first frame is frame number 0 so that end is M-1, e.g. for 10 frames, end=9.
i/f
 This is the option which determines in which domain the fidelity constraint is computed. This is a binary flag with 1 for the image domain and 0 for the Fourier domain. The advantage of the image domain is that a bad pixel mask can be used, which is input as the data support conv_sup. The advantage of the Fourier domain is that the convergence can be stabilized by using a weighting on the different spatial frequencies which can prevent noise amplification. 

  (Note: Bad pixel rejection via conv_sup is NOT supported in Fourier domain.) 

conv_wf
  This is a 2-D [or 3_D] Wiener Filter used as a weighting function in conv-spectrum space. 

  If the corresponding Wiener filter is not available, you must provide a default file of size N x N (conv_wf) = 1.0.

sky
 The algorithm always expects a sky file even when the sky has been subtracted from the data. The sky can vary from frame-to-frame so the data array has a size of N x N x M. If the sky level is constant across all frames, then a single frame, i.e. N x N, sky file can be used. For sky subtracted data, the array should be set to zero.  Note that the sky image is not modifiedby the code, but a multiplicative factor (Eta) is. This value is minimized by the routine just like the other object and psf variables. The code trys to fit that sky frame to the DC level of the convolution image(s) in a least-squares sense.

Convergence

How do you know when you've converged? This is a difficult question for many of these non-linear iterative algorithms. It's possible to overfit the data resulting in noise amplification which can significantly affect the reconstruction although still minimizing the combined error-metric. To address this issue, we have implemented a scheme which will terminate the convergence when the squared residual of the measured and estimated data reached the noise level.

For most data, the noise can be modeled by a zero-mean Gaussian characterized by a specific rms value, i.e standard deviation. This can be estimated from empty regions from the data frames or from corresponding sky or dark frames. This noise value needs to be added to the FITS header of the convolution file so that the algorithm knows when to stop. The keyword for this value is N0. This assumes that all of the frames have the same SNR which is typically realistic in that frames with the same exposure time will have the same noise statistics and also the same total power.

Computing The Noise Limit

Remember that the input data is normalized on a frame-by-frame basis such that its total power or volume is unity. Thus the measured rms noise also has to normalized by the same factor. In the current code set-up this is calculated internally. Thus the user need only enter N0.

Example

A sample data set has been created and can be downladed from here. It comprises 10 frames (64 x 64 pixels) of speckle interferometric observations of a multiple star system. It is advised that the first-time idac user test out the algorithm on this data set and compare the results to those found here. This page illustrates the application of idac to this sample data.