SpeX Iterated PSFs
Sky Coyote, November 2004

We have as input data a movie (data0002.movie) of 8192 64x64 10ms exposure images of a star taken with the IRTF SpeX guider by Doug Toomey on 5/27/00. See the log file for more details. These images show speckle interference of the star passed through the atmosphere and the PSF of the telescope. A typical frame of this data is shown below:

The goal is to use these images to construct a good PSF of the telescope. The procedure is to iterate on a PSF, using as a starting point the shift-and-add best images from the movie.

We assume the following:

   I(t) = (O * A(t)) * P
        = A(t) * P
where I(t) = data image frame of movie
         O = object, a point source
      A(t) = atmosphere PSF of each frame
         P = constant PSF of telescope
         * = convolution

This convolution process is illustrated below, where the boxes outlined in grey are known quantities, and P is to be determined:

The iterated procedure is:

   A(t) = I(t) *' P
   P(t) = I(t) *' A(t)
     P' = f({P(t)})
where P(t) = 'partial' telescope PSF of each frame
        P' = new PSF
         f = some function of all P(t)
        *' = deconvolution via Lucy-Richardson

As a comparison/alternative to this procedure, we also use the autocorrelation procedure from Labeyrie 1970 [1]:

                    I(t) = O * (abs(F(PP)))^2
                         = (abs(F(PP)))^2

                 F(I(t)) = AC(PP)
        (abs(F(I(t))))^2 = (abs(AC(PP)))^2
   Sum((abs(F(I(t))))^2) = Sum((abs(AC(PP)))^2)
where PP = telescope pupil
       F = Fourier transform
     abs = absolute value (magnitude)
      AC = autocorrelation
     Sum = summation

The overall procedure is contained in the IDL program psfiter.pro. The steps are:

1. A 256 frame subset of the 8192 frame movie was extracted and saved in data.movie. A 'best' PSF was created by registering all 8192 frames of the movie so that the brightest pixel was in the center of the frame, sorting the registered frames based on the value of the brightest pixel, and then summing the top 500 frames. A starting PSF was created by shifting and adding all 256 frames of the subset movie. These PSFs, normalized to a total of 1, are presented below:

The goal of this test was to see if the iterated procedure would transform the starting PSF into something that looked like the 'best' PSF.

2. 10 iterations of the procedure were used to create a series of PSFs. At each step, 10 iterations of Lucy-Richardson deconvolution were performed. The function f({P(t)}) is the simple sum of all P(t). The first few results, normalized to a total of 1, are shown below (steps 0-5):

However, the iterated PSFs do not approach the 'best' PSF, but instead appear to approach a delta distribution. This is an unfortunate consequence of the fact that a delta distribution is a trivial solution of the equation I(t) = A(t) * P, in which case the A(t) should also approach the input data I(t). This is indeed the case, as can be seen by a comparison of one frame of I(t) and the final A(t):

3. The Labeyrie autocorrelation was calculated and compared to the final value of P. In this case, there is good agreement between the two results:

also shown here at a stretched scale:

Unfortunately, the Labayrie autocorrelation does not include phase information, so that it is not possible to reconstruct a non-symmetric PSF from this pupil function. Nevertheless, we believe that placing an appropriate constraint on the function f({P(t)}), or on the intermediate values A(t), will prevent the solution P from degrading to the trivial result.

We have found that other methods of constructing a sequence of better PSFs (e.g. evolution of a population of PSFs via a genetic algorithm) also converge toward a delta distribution without the use of constraints. Possible constraints on A(t) include that they be non-negative and real. However, in the case of the above test, all iterated values of A(t) were non-negative and real, due to the Lucy-Richarson deconvolution method used. Obviously, further investigation and experimentation is required.

[1] 'Attainment of Diffraction Limited Resolution in Large Telescopes by Fourier Analysing Speckle Patterns in Star Images', A. Labeyrie, Astron. & Astrophys. 6, 85-87 (1970).