2D PSF from 1D component

13 次查看（过去 30 天）
SvB
SvB 2018-2-28
编辑： Iryna Fizor 2019-11-26
This is something that has been annoying me to no end today, and I'm not sure whether it's because of a lack of understanding the math, not being able to translate my idea to matlab or that everything actually does work and it just doesn't get better than this. I hope someone can help me out!
Long story short: Say I have a 2D MTF (Modulation Transfer Function) that I transform into a 2D PSF (Point Spread Function) using the ifft2() function in Matlab. This works fine, but I'll need to do many of these calculations at high resolution (say, 4000*6000), so using ifft2() will take up a lot of time. I figured that if the MTF / PSF were circularly symmetrical, I should be able to do a 1D transformation from, say, MTF_x into PSF_x, and then reconstruct the full 2D PSF from PSF_x. However, it just doesn't seem to work (right).
Please note that while the example here uses a normal distribution, the actual MTF / PSFs are not. (They are circularly symmetric, though). Therefore, it won't be possible to extract parameters such as center and sigma and use them to plot a new Gauss as PSF.
A working example code for the 2D approach:
% Preparations
mu = [0]; % Center of MTF
Sigma = [2 2;]; % Sigma along X and Y axis
fx = -100:1:100; % Frequency coordinates X-axis
fy = -100:1:100; % Frequency coordinates Y-axis
[X,Y] = meshgrid(fx,fy); % Full mesh-grid
% Calculating MTF
MTF2 = mvnpdf([X(:) Y(:)],mu,Sigma); % Normal distribution along X and Y axis
MTF2 = reshape(MTF2,length(fy),length(fx)); % Reshape into a matrix.
% Calculating PSF
PSF2 = abs(fftshift(ifft2(fftshift(MTF2)))); % Calculate PSF
PSF2 = PSF2./sum(PSF2(:)); % Normalise PSF
This yields a nice 2D PSF. Now, taking the 1D approach:
% Calculating MTF (1D)
MTF1=MTF2(:,ceil(length(fx)/2)); % Take the centerline of the 2D MTF to be the 1D MTF
% Alternatively, you can create the MTF rather than taking the centerline by using:
% MTF1 = normpdf(fx,0,sqrt(2)); % Center 0, sigma is sqrt(2), I checked and it fits the approach taken above.
% Calculating PSF (1D)
PSF1=abs(fftshift(ifft(fftshift(MTF1))));
% Conversion from 1D to circularly symmetric 2D
r = sqrt(X.^2+Y.^2); % Generate matrix with radial coordinates
PSF1_R = interp1(fx,PSF1,r,'cubic'); % For each value of r in the matrix, find the corresponding PSF1-value.
PSF1_R = PSF1_R./nansum(PSF1_R(:)); % Normalise. NaNs are present (and ignored) because r is larger than X in the corners of the image.
Then, I plot the ratio of these two results as follows:
imshow(PSF2./PSF1_R,[0 2])
Ideally, this should be a perfect grey image since all values should be 1. This situation is not ideal, since r>x and thus interpolation is not possible in the 4 edges of the image, so you'd expect a perfect grey circle instead. At first glance it doesn't look too bad, however, focussing more on the values near 1 yields:
imshow(PSF2./PSF1_R,[0.98 1.02])
Which certainly isn't a nice flat circle.
Question 1) Are these rounding errors, or what is causing this? I'd say that by calculating r the way I did, and interpolating the PSF for r the way I did, I should be getting a circularly symmetric shape, yet I'm not.
Now, as a bonus: When defining fx and fy, reduce the step size for both to 0.1 (from 1). I would expect the results to improve, since you work within the same (frequency) range, but at smaller step sizes and thus higher spatial resolution in the PSF. However, instead the result is something that wouldn't look bad as an easter egg.
Question 2) What's going on here? What is causing this?
Note A) I just realise that while fx and fy are in the frequency domain (matched to the MTF), r is in the space domain (matched to PSF). However, I do use interp1() on fx and r together. Might this be causing issues?
Note B) Using linear rather than cubic interpolation in interp1() will generate an additional 'structure' in the shown images, which leads me to think that Question 1) should at least partially be caused by rounding errors, but it doesn't explain the lobes.
  1 个评论
Iryna Fizor
Iryna Fizor 2019-11-26
编辑：Iryna Fizor 2019-11-26
Hello SvB, your research on this is very interesting. I will be very appreciated if you have found the answer or hints to your questions and could share your share them. I need to become an approximate 2d PSF from 1D MTF.

请先登录，再进行评论。

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

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Uniform Distribution (Continuous) 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by