Calculate encircled energy from point spread function

Hi all,

I have a point spread function and I am trying to reproduce the encircled energy (EE) plot to find the encircled energy of it (see example of what I would like). I have attempted to calculate the EE and plot the EE against the radial distance from the centre of the PSF. However, the EE should be from 0 to 1 on the y axis but in my plot (see attached), it doesn’t seem to be the case. I have attached the code that I have used:

horizontalProfile = psf(128, :);
cdf = cumsum(horizontalProfile(128:end));
% Normalize
cdf = cdf / sum(cdf);
plot(cdf)

The value 128 is the location of the centre of the PSF (grid size 256 x 256). If anyone could kindly look at it and tell me what I am doing wrong, I would be so grateful! I have also attached the image of the PSF generated.

[Vanessa]

Hi David,

Thank you so much for your answer, the graph is now going from 0 to 1. However, it starts at 0.2, is this okay?

[David Goodmanson]

Hi Vanessa, I withdrew the answer. It was not quite right because you have 2-d situation. And that is related to why you are getting a start of .2.

The 256x256 psf matrix, do the rows and columns correspond to variables r and theta, or to variables x and y? I would guess it's x and y, in which case it takes some work to get the answer unless you use the approximation that the psf function is axially symmetric and does not depend on theta. I will repost an answer in awhile.

[Vanessa]

Hi David,

Thank you so much for your help, it is much appreciated and you guessed right! the rows and column correspond to x and y.

Answer 1

[David Goodmanson]

Hi Vanessa,

I will assume that the psf is a function of intensity and not amplitude; otherwise you would have to square the psf in order to get intensity. In the posted picture the horizontal axis is in minutes, implying something is changing with time, but let's assume that the intensity is uniform in time. This means using energy/sec rather than energy but does not really change the cdf calculation, which is a spacial calculation.

Energy/sec passing through the surface is intensity times times area. The psf is axially symmetric, so for an increment dr, the area da is simply the annulus 2*pi*r*dr. The dr part is associated with the array spacing in the psf array and the 2*pi will get taken care of by normalization. That leaves r, so to create the pdf you need to multiply the psf by r. Since the array has an even number of points it seems like the center could be either 128 or 129 (or 128.5). But let's assume 128 as you did.

It's the pdf that gets normalized to 1, not the cdf.

The following code uses the same data from the psf as you did.

horizontalProfile = psf(128, :);
HP = horizontalProfile(128:end);
r = 0:length(HP)-1;            % 129 points
HPpdf = HP.*r;
HPpdf = HPpdf/sum(HPpdf);      % normalize
cdf = cumsum(HPpdf);
plot(r,cdf)

Here the dr spacing was taken as 1, so you might have to rescale the r axis if you know how the psf was sampled.

[Vanessa]

HI David,

Thank you for your code, I tested it out with:

psf = fspecial('gaussian', 48, 5) ;

and got the right graph with it: