karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Working version of the influence of aberrations on a PSF

Browse Source
This commit is contained in:
Karthik 2025-02-23 01:24:32 +01:00
parent 66576da594
commit 612bae7fa1
1 changed files with 62 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

128
Estimations/ModellingImagingAberrations.m
View File

@ -13,24 +13,27 @@ plotZernike(3, -1, NumberOfGridPoints) % Coma (Z3^-1, y-direction)
plotZernike(4, 0, NumberOfGridPoints) % Spherical Aberration (Z4^0)
%% Compute and plot Aberrated PSF, Image
NumberOfGridPoints = 1024; % Number of grid points per side
PupilRadius = 0.010; % Radius of pupil [m]
Length = 0.5; % Total size of the grid [m]
NumberOfGridPoints = 256; % Number of grid points per side
PupilRadius = 10E-3; % Radius of pupil [m]
Length = 32E-6; % Total size of the grid [m]
GridSpacing = Length / NumberOfGridPoints; % Grid spacing [m]
Wavelength = 421e-9; % Optical wavelength [m]
ImageDistance = 0.7; % Image distance [m]
Wavelength = 421E-9; % Optical wavelength [m]
FocalLength = 3*PupilRadius; % Focal Length
% Generate PSF
C = [0.3, 0.2, 0.1, -0.1, 0.4]; % Zernike coefficients
[X, Y, PSF] = generateAberratedPSF(C, PupilRadius, NumberOfGridPoints, GridSpacing);
plotAberratedPSF(X, Y, PSF);
C = [1.5, 0.0, 0.1, -0.1, 0.5]; % Zernike coefficients
PSF = generateAberratedPSF(C, Wavelength, PupilRadius, NumberOfGridPoints, GridSpacing, FocalLength);
xvals = (-NumberOfGridPoints/2 : NumberOfGridPoints/2-1) * GridSpacing * 1E6;
yvals = (-NumberOfGridPoints/2 : NumberOfGridPoints/2-1) * GridSpacing * 1E6;
plotImage(xvals, yvals, PSF, 'PSF');
% Generate object
Object = generateObject(Wavelength, ImageDistance, NumberOfGridPoints, GridSpacing);
% Object = generateObject(Wavelength, NumberOfGridPoints, GridSpacing, FocalLength);
% Convolve the object with the PSF to simulate imaging
Image = convolveObjectWithPSF(abs(Object).^2, PSF, 1);
% Image = convolveObjectWithPSF(abs(Object).^2, PSF);
% plotImage(xvals, yvals, Image, 'Image of Object formed by convolving with the PSF');
%% Functions
function Z = computeZernikePolynomials(n, m, r, theta)
@ -46,19 +49,19 @@ function Z = computeZernikePolynomials(n, m, r, theta)
if n == 2 && m == 0
% Defocus (Z2^0)
Z = 2 * r.^2 - 1;
Z = sqrt(3) * (2 * r.^2 - 1);
elseif n == 2 && m == 2
% Astigmatism (Z2^2)
Z = r.^2 .* cos(2 * theta);
% Vertical Astigmatism (Z2^2)
Z = sqrt(6) * r.^2 .* cos(2 * theta);
elseif n == 3 && m == 1
% Coma (Z3^1)
Z = (3 * r.^3 - 2 * r) .* cos(theta);
% Horizontal Coma (Z3^1)
Z = sqrt(8) * (3 * r.^3 - 2 * r) .* cos(theta);
elseif n == 3 && m == -1
% Coma (Z3^-1)
Z = (3 * r.^3 - 2 * r) .* sin(theta);
% Vertical Coma (Z3^-1)
Z = sqrt(8) * (3 * r.^3 - 2 * r) .* sin(theta);
elseif n == 4 && m == 0
% Spherical Aberration (Z4^0)
Z = 6 * r.^4 - 6 * r.^2 + 1;
Z = sqrt(5) * (6 * r.^4 - 6 * r.^2 + 1);
else
% Default to zero if no known Zernike polynomial matches
Z = 0;
@ -94,7 +97,7 @@ function plotZernike(n, m, NumberOfGridPoints)
shading interp;
end
function [X, Y, PSF] = generateAberratedPSF(C, PupilRadius, NumberOfGridPoints, GridSpacing)
function PSF = generateAberratedPSF(C, Wavelength, PupilRadius, NumberOfGridPoints, GridSpacing, FocalLength)
% C is the vector of Zernike coefficients [C_defocus, C_astigmatism, C_coma, C_spherical, ...]
% NumberOfGridPoints is the number of points for the grid (NxN grid)
% PupilRadius is the radius of the pupil aperture
@ -103,72 +106,65 @@ function [X, Y, PSF] = generateAberratedPSF(C, PupilRadius, NumberOfGridPoints,
[X, Y] = meshgrid((-NumberOfGridPoints/2 : NumberOfGridPoints/2-1) * GridSpacing);
% Convert (x, y) to polar coordinates (r, theta)
[theta, r] = cart2pol(X, Y);
[theta, ~] = cart2pol(X, Y);
% Pupil function: 1 inside the pupil radius, 0 outside
P = generateCircularMask(X, Y, 2 * PupilRadius); % 2D mask
% Pupil function
P = generatePupil(Wavelength, PupilRadius, NumberOfGridPoints, GridSpacing, FocalLength);
% Wavefront error from Zernike polynomials
W = C(1) * computeZernikePolynomials(2, 0, r, theta) + ... % Defocus (Z2^0)
C(2) * computeZernikePolynomials(2, 2, r, theta) + ... % Astigmatism (Z2^2)
C(3) * computeZernikePolynomials(3, 1, r, theta) + ... % Coma (Z3^1, x-direction)
C(4) * computeZernikePolynomials(3, -1, r, theta) + ... % Coma (Z3^-1, y-direction)
C(5) * computeZernikePolynomials(4, 0, r, theta); % Spherical Aberration (Z4^0)
W = C(1) * computeZernikePolynomials(2, 0, P, theta) + ... % Defocus (Z2^0)
C(2) * computeZernikePolynomials(2, 2, P, theta) + ... % Astigmatism (Z2^2)
C(3) * computeZernikePolynomials(3, 1, P, theta) + ... % Coma (Z3^1, x-direction)
C(4) * computeZernikePolynomials(3, -1, P, theta) + ... % Coma (Z3^-1, y-direction)
C(5) * computeZernikePolynomials(4, 0, P, theta); % Spherical Aberration (Z4^0)
W(P>1) = 0;
E = exp(1i*2*pi*W); % Complex amplitude
E(P>1) = 0; % Impose aperture size
% Fourier transform of the pupil function with aberrations
PSF = abs(calculateExtendedFFT2(P .* exp(-1i * 2*pi * W), GridSpacing)).^2;
PSF = abs(calculateFFT2(E)).^2;
PSF = PSF/sum(PSF(:));
end
function plotAberratedPSF(X, Y, PSF)
% C: Zernike coefficients [C_defocus, C_astigmatism, C_coma_x, C_coma_y, C_spherical]
% PupilRadius: Radius of the pupil aperture
% NumberOfGridPoints: Number of points for plotting
% Generate PSF using the updated modelPSF function
PSF = PSF / max(PSF(:));
figure(1)
function plotImage(xvals, yvals, image, titlestring)
figure
clf
set(gcf,'Position',[50 50 950 750])
surf(X, Y, PSF, 'EdgeColor', 'none');
xlim([-0.05 0.05])
ylim([-0.05 0.05])
view(2); % 2D view
shading interp;
imagesc(xvals, yvals, image);
colorbar;
colormap jet;
title('PSF', 'FontSize', 16);
title(titlestring, 'FontSize', 16);
xlabel('X', 'FontSize', 16);
ylabel('Y', 'FontSize', 16);
end
function C = convolveObjectWithPSF(A, B, GridSpacing)
N = size(A, 1);
C = calculateExtendedInverseFFT2(calculateExtendedFFT2(A, GridSpacing) .* calculateExtendedFFT2(B, GridSpacing), 1/(N*GridSpacing));
function P_Norm = generatePupil(Wavelength, PupilRadius, NumberOfGridPoints, GridSpacing, FocalLength)
[X,Y] = meshgrid((-NumberOfGridPoints/2 : NumberOfGridPoints/2-1) * ((Wavelength * FocalLength) / (NumberOfGridPoints * GridSpacing)));
P = sqrt(X.^2 + Y.^2);
P_Norm = P/PupilRadius; % Pupil normalized by aperture radius
assert(max(P_Norm(:))>=sqrt(2),'Sampling is not sufficient to reconstruct the entire wavefront.');
end
function z = generateCircularMask(x, y, D)
r = sqrt(x.^2+y.^2);
z = double(r<D/2);
z(r==D/2) = 0.5;
function G = calculateFFT2(g)
G = fftshift(fft2(ifftshift(g)));
end
function G = calculateExtendedFFT2(g, Delta)
G = fftshift(fft2(fftshift(g))) * Delta^2;
%{
function g = calculateInverseFFT2(G)
g = ifftshift(ifft2(fftshift(G)));
end
function g = calculateExtendedInverseFFT2(G, Delta_f)
N = size(G, 1);
g = ifftshift(ifft2(ifftshift(G))) * (N * Delta_f)^2;
end
function ret = generateRectangle(x, D)
if nargin == 1, D = 1; end
x = abs(x);
ret = double(x<D/2);
ret(x == D/2) = 0.5;
end
function obj = generateObject(Wavelength, ImageDistance, NumberOfGridPoints, GridSpacing)
function obj = generateObject(Wavelength, NumberOfGridPoints, GridSpacing, FocalLength)
% Image-plane coordinates
[U, V] = meshgrid((-NumberOfGridPoints/2 : NumberOfGridPoints/2-1) * ((Wavelength * ImageDistance) / (NumberOfGridPoints * GridSpacing)));
obj = (generateRectangle((U-1.4e-4)/5e-5) + generateRectangle(U/5e-5)+ generateRectangle((U+1.4e-4)/5e-5)) .* generateRectangle(V/2e-4);
[U, V] = meshgrid((-NumberOfGridPoints/2 : NumberOfGridPoints/2-1) * ((Wavelength * FocalLength) / (NumberOfGridPoints * GridSpacing)));
obj = ;
end
function C = convolveObjectWithPSF(A, B)
C = calculateInverseFFT2(calculateFFT2(A) .* calculateFFT2(B));
end
%}