Calculate encircled energy from point spread function

9 次查看（过去 30 天）
Vanessa
Vanessa 2017-9-9
编辑： Vanessa 2017-9-11
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.
  3 个评论
显示 1更早的评论
David Goodmanson
David Goodmanson 2017-9-9
编辑：David Goodmanson 2017-9-9
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
Vanessa 2017-9-10
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

David Goodmanson
David Goodmanson 2017-9-10
编辑：David Goodmanson 2017-9-10
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.
  1 个评论
Vanessa
Vanessa 2017-9-11
编辑：Vanessa 2017-9-11
HI David,
Thank you for your code, I tested it out with:
psf = fspecial('gaussian', 48, 5) ;
and got the right graph with it:

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Lighting, Transparency, and Shading 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by