diff --git a/Estimations/ModellingImagingAberrations.m b/Estimations/ModellingImagingAberrations.m index b5bef97..98ddd95 100644 --- a/Estimations/ModellingImagingAberrations.m +++ b/Estimations/ModellingImagingAberrations.m @@ -1,22 +1,37 @@ -%% Zernike Polynomials +%% Compute and plot Zernike Polynomials -resolution = 100; % Resolution +NumberOfGridPoints = 100; % Resolution -plotZernike(2, 0, resolution) % Defocus (Z2^0) +plotZernike(2, 0, NumberOfGridPoints) % Defocus (Z2^0) %% -plotZernike(2, 2, resolution) % Astigmatism (Z2^2) +plotZernike(2, 2, NumberOfGridPoints) % Astigmatism (Z2^2) %% -plotZernike(3, 1, resolution) % Coma (Z3^1, x-direction) +plotZernike(3, 1, NumberOfGridPoints) % Coma (Z3^1, x-direction) %% -plotZernike(3, -1, resolution) % Coma (Z3^-1, y-direction) +plotZernike(3, -1, NumberOfGridPoints) % Coma (Z3^-1, y-direction) %% -plotZernike(4, 0, resolution) % Spherical Aberration (Z4^0) +plotZernike(4, 0, NumberOfGridPoints) % Spherical Aberration (Z4^0) -%% Aberrated PSF -C = [0.0, 0.0, 0.0, 0.0, 0.0]; % Zernike coefficients -pupil_radius = 0.5; % Pupil radius (normalized) +%% 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] +GridSpacing = Length / NumberOfGridPoints; % Grid spacing [m] +Wavelength = 421e-9; % Optical wavelength [m] +ImageDistance = 0.7; % Image distance [m] -plotAberratedPSF(C, pupil_radius, resolution); +% 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); + +% Generate object +Object = generateObject(Wavelength, ImageDistance, NumberOfGridPoints, GridSpacing); + +% Convolve the object with the PSF to simulate imaging +Image = convolveObjectWithPSF(abs(Object).^2, PSF, 1); + +%% Functions function Z = computeZernikePolynomials(n, m, r, theta) % Zernike polynomial function for radial and angular coordinates (r, theta) @@ -50,13 +65,13 @@ function Z = computeZernikePolynomials(n, m, r, theta) end end -function plotZernike(n, m, resolution) +function plotZernike(n, m, NumberOfGridPoints) % n: radial order % m: azimuthal frequency - % resolution: number of points for plotting + % NumberOfGridPoints: number of points for plotting % Create a grid of (r, theta) coordinates - [theta, r] = meshgrid(linspace(0, 2*pi, resolution), linspace(0, 1, resolution)); + [theta, r] = meshgrid(linspace(0, 2*pi, NumberOfGridPoints), linspace(0, 1, NumberOfGridPoints)); % Calculate the Zernike polynomial for the given (n, m) Z = computeZernikePolynomials(n, m, r, theta); @@ -79,16 +94,19 @@ function plotZernike(n, m, resolution) shading interp; end -function [theta, r, PSF] = modelPSF(C, pupil_radius, resolution) +function [X, Y, PSF] = generateAberratedPSF(C, PupilRadius, NumberOfGridPoints, GridSpacing) % C is the vector of Zernike coefficients [C_defocus, C_astigmatism, C_coma, C_spherical, ...] - % resolution is the number of points for the grid (NxN grid) - % pupil_radius is the radius of the pupil aperture + % NumberOfGridPoints is the number of points for the grid (NxN grid) + % PupilRadius is the radius of the pupil aperture - % Create a grid of (r, theta) coordinates - [theta, r] = meshgrid(linspace(0, 2*pi, resolution), linspace(0, 1, resolution)); + % Create a grid of (x, y) of pupil-plane coordinates in Cartesian space + [X, Y] = meshgrid((-NumberOfGridPoints/2 : NumberOfGridPoints/2-1) * GridSpacing); + % Convert (x, y) to polar coordinates (r, theta) + [theta, r] = cart2pol(X, Y); + % Pupil function: 1 inside the pupil radius, 0 outside - P = double(r <= pupil_radius); % 2D mask + P = generateCircularMask(X, Y, 2 * PupilRadius); % 2D mask % Wavefront error from Zernike polynomials W = C(1) * computeZernikePolynomials(2, 0, r, theta) + ... % Defocus (Z2^0) @@ -98,26 +116,22 @@ function [theta, r, PSF] = modelPSF(C, pupil_radius, resolution) C(5) * computeZernikePolynomials(4, 0, r, theta); % Spherical Aberration (Z4^0) % Fourier transform of the pupil function with aberrations - PSF = abs(fftshift(fft2(P .* exp(1i * W)))).^2; % Intensity distribution + PSF = abs(calculateExtendedFFT2(P .* exp(-1i * 2*pi * W), GridSpacing)).^2; end -function plotAberratedPSF(C, pupil_radius, resolution) +function plotAberratedPSF(X, Y, PSF) % C: Zernike coefficients [C_defocus, C_astigmatism, C_coma_x, C_coma_y, C_spherical] - % pupil_radius: Radius of the pupil aperture - % resolution: Number of points for plotting - + % PupilRadius: Radius of the pupil aperture + % NumberOfGridPoints: Number of points for plotting % Generate PSF using the updated modelPSF function - [theta, r, PSF] = modelPSF(C, pupil_radius, resolution); - - % Convert polar to Cartesian coordinates for plotting - [X, Y] = pol2cart(theta, r); - + PSF = PSF / max(PSF(:)); figure(1) 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 - axis equal tight; shading interp; colorbar; colormap jet; @@ -125,3 +139,36 @@ function plotAberratedPSF(C, pupil_radius, resolution) 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)); +end + +function z = generateCircularMask(x, y, D) + r = sqrt(x.^2+y.^2); + z = double(r