The image below demonstrates that a better PSF can be found using a single step of the iterative process:
The windows show the following:
Top left (0): 1 frame of data, normalized to 1 Middle left (1): Initial PSF, normalized to 1 Bottom (2): a = lucy(data, psf, n=50), normalized to 1 Middle center (3): p = lucy(data, a, n=50), normalized to 1 Middle right (4): Point source = 1 Top 2nd (5): d2 = convolve(a, psf) Top 3rd (6): d3 = convolve(a, p) Top right (7): d4 = convolve(a, point)3 residuals were calculated as:
r2 = total((data - d2)^2) = 0.000320490 r3 = total((data - d3)^2) = 5.85313e-05 r4 = total((data - d4)^2) = 0.00611659r3 is the least of these values, indicating that p is a better PSF than either the initial PSF or a point source. Therefore, there exists a 'best' PSF somewhere between the initial PSF and a delta distribution. Assuming that this result holds for all frames of input data, the question is then how to combine the p(t) to produce a composite PSF which is as good as, or better than, all of them.