Effacer les filtres
Effacer les filtres

Error in Generating Point Spread Function

18 vues (au cours des 30 derniers jours)
Asser Abdelgawad
Asser Abdelgawad le 6 Juin 2022
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?

Réponses (1)

Debadipto
Debadipto le 6 Fév 2024
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:

Produits


Version

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by