The code I was written by Jason Marshall, he wrote an SPIE paper that
describes the algorithm in part (see
http://ao.jpl.nasa.gov/papers/scientific/4007-731.pdf)

Estimate of Strehl:
We start by assuming the images have been bad pixel corrected, background
subtracted and flat fielded. We then use a C program called SExtractor,
which as I understand it is a photometric code, designed to work in crowded
fields. This program provides an estimate of the image background, total
flux from the selected object and refined x,y location.

We then generate a airy function centered on the same sub-pixel location as
the science image at the same pixel scale (arc/pixel) as the science image.
The airy function generation includes any secondary obscuration, but not
spiders. The airy function is then scaled to have the same total flux as the
science image.

At this point we have a science image and "perfect" image, with the same
total flux, on the same plate scale and centered on the same sub-pixel.

Next we calculate the peak intensity of each image using IDL's cubic
interpolate function, with a factor of -0.5. The Strehl is just the ratio of
the two peak intensities.

Random and systematic error of Strehl estimate:

We generate a "pseudo" AO image, which has the same Strehl and peak
intensity as the measured values. To do this the user must enter the seeing.
This is essentially a diffraction limited image on top of a halo that is
lambda/r0. This "pseudo" image is then inserted into a empty portion of the
original science image and then Strehl calculated as above. This is done for
5 random locations. The difference between the mean of these 5 Strehls and
the estimated Strehl is called the Systematic error. The RMS of the 5
Strehls is called the random error.