Here's how my Strehl meter works.

The data in the neighborhood of the center of the image is extracted. The
intensity of the image is oversampled by a factor of 8 (but this can be
changed) using FFT interpolation. Then the maximum of the value of this
oversampled data is found. The maximum value is then divided by the total
intensity over a circular aperture (typically of 1" radius). The
normalization factor needed to define the Strehl is found by running the
diffraction-limited image through the Strehl meter and setting the Strehl
to be unity.

The diffraction-limited PSF is found as follows.
Using discrete Fourier transforms, the perfect PSF is calculated on a four
times oversampled grid and then rebinned to the right sampling. In this
way, sampling vs binning issues are minimized.

I would like to point out that the code is completely insensitive to
subpixel shifts for Nyquist sampled data (this is why I wrote the
Strehl meter in the first place). It works extremely well for
high Strehls when Nyquist sampled. It understimates the Strehls
when the Strehls are low as the true maximum cannot be found using
interpolation due to the loss of information implicit in the
pixellation. Maybe a solution to this problem would be to use a
pre-computed "fudge factor" as a function of the measured Strehl to boost
the low Strehls.