diff --git a/Estimations/ModellingImagingAberrations.m b/Estimations/ModellingImagingAberrations.m index 98ddd95..aac55e1 100644 --- a/Estimations/ModellingImagingAberrations.m +++ b/Estimations/ModellingImagingAberrations.m @@ -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