Error in Generating Point Spread Function

조회 수: 6 (최근 30일)
Asser Abdelgawad
Asser Abdelgawad 2022년 6월 6일
답변: Debadipto 2024년 2월 6일
After finding the Modulation Transfer Function, I used this method to rotate it about its vertical axis to obtain the 2D MTF. I then performed ifft2() on that to generate the PSF. However, my PSF looks off for some reason.
%restore fitted MTF properties
freq = q/T;
MTF_fit = abs(fft(yLSF_fit));
MTF_fit = fftshift(MTF_fit);
f = fit(freq.', MTF_fit.', 'gauss2');
% rotate MTF in fourier domain
n = length(freq);
x0 = freq;
y0 = f(x0).';
theta = linspace(0,2*pi,n).';
x = x0.*cos(theta);
y = x0.*sin(theta);
z = repmat(y0,[n,1]);
MTF2 = z;
%obtain 2D PSF using ifft2 on 2D MTF
PSF2 = abs(ifftshift(ifft2((MTF2)); % Calculate PSF
% or PSF2 = abs(fftshift(ifft2(fftshift(MTF2))));
PSF2_normalized = (PSF2-min(PSF2(:)))./(max(PSF2(:))-min(PSF2(:))); % Normalize PSF
%plot 2D MTF, PSF, and normalized PSF
subplot(3,2,4)
surf(x,y,z)
xlabel pixels
ylabel pixels
zlabel amplitude
title('2D MTF')
shading interp
axis tight
subplot(3,2,5)
surf(PSF2);
title('PSF')
shading interp
axis tight
xlabel pixels
ylabel pixels
zlabel amplitude
subplot(3,2,6)
surf(PSF2_normalized);
title('PSF (Normalized)')
shading interp
axis tight
rotate3d
xlabel pixels
ylabel pixels
zlabel amplitude
rotate3d('on');
What could be going wrong? I suspect it has something to do with the method I used for rotation because it generates exact copies of the same number that repeats throughout each column of the matrix, where each column holds a different value. See the MTF2 variable in the .mat workspace I have attached.
I've also attached another workspace (validMTF2.mat) with a correct MTF2 (i.e: one that has a viable shape and actually generates a 3D Gaussian PSF when ifft2() is applied) for reference. In this matrix, the columns aren't comprised of a single value, which is why I am suspicious of the rotation method I used, even though its graph looks perfect. Beware the values are very small though so they might show up as zeros. Thanks for the help!
Edit: here is what happens when n = 10 in the above code. The PSF has a more 3D shape in this case (higher n values make it thinner and thinner, which I don't understand). The tradeoff here is a less accurate 2D MTF shape as shown by the spikiness. How can I fix this issue and maintaing a high n value for accuracy?

답변 (1개)

Debadipto
Debadipto 2024년 2월 6일
Replicating the 1D MTF using repmat actually creates a 2D MTF with repeated values, which is not rotationally symmetric. Instead, you should create a 2D grid and evaluate the 1D MTF at each point's radial distance from the center.
The following approach should give you a rotationally symmetric 2D MTF:
%restore fitted MTF properties
freq = q/T
MTF_fit = abs(fft(yLSF_fit));
MTF_fit = fftshift(MTF_fit);
f = fit(freq.', MTF_fit.', 'gauss2');
% Define the grid size
n_val = 100;
% Create a 2D grid of points
[x, y] = meshgrid(linspace(-max(freq), max(freq), n_val), linspace(-max(freq), max(freq), n_val));
% Calculate the radial distance from the center for each point
r = sqrt(x.^2 + y.^2);
% Evaluate the 1D MTF at each radial distance to get a 2D rotationally symmetric MTF
% Using arrayfun to apply 'f' to each element of 'r'
z = arrayfun(@(r_val) f(r_val), r);
MTF2 = z;
%obtain 2D PSF using ifft2 on 2D MTF
PSF2 = abs(ifftshift(ifft2((MTF2)))); % Calculate PSF
% or PSF2 = abs(fftshift(ifft2(fftshift(MTF2))));
PSF2_normalized = (PSF2-min(PSF2(:)))./(max(PSF2(:))-min(PSF2(:))); % Normalize PSF
%plot 2D MTF, PSF, and normalized PSF
subplot(3,2,4)
surf(x,y,z)
xlabel pixels
ylabel pixels
zlabel amplitude
title('2D MTF')
shading interp
axis tight
subplot(3,2,5)
surf(PSF2);
title('PSF')
shading interp
axis tight
xlabel pixels
ylabel pixels
zlabel amplitude
subplot(3,2,6)
surf(PSF2_normalized);
title('PSF (Normalized)')
shading interp
axis tight
rotate3d
xlabel pixels
ylabel pixels
zlabel amplitude
rotate3d('on');
Here is the plot generated on running the code:

카테고리

Help CenterFile Exchange에서 Surface and Mesh Plots에 대해 자세히 알아보기

제품


릴리스

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by