diff --git a/Data-Analyzer/plotImages.m b/Data-Analyzer/plotImages.m index f7d01e1..4675378 100644 --- a/Data-Analyzer/plotImages.m +++ b/Data-Analyzer/plotImages.m @@ -1,12 +1,12 @@ %% Parameters -groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", "/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", "/images/Vertical_Axis_Camera/in_situ_absorption"]; +groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", "/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", "/images/Vertical_Axis_Camera/in_situ_absorption"]; -folderPath = "C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/20/"; +%{ -run = '0140'; +folderPath = "C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/"; -folderPath = strcat(folderPath, run); +run = '0140'; cam = 4; @@ -16,7 +16,25 @@ span = [50, 50]; fraction = [0.1, 0.1]; pixel_size = 5.86e-6; -removeFringes = true; +removeFringes = false; + +%} + +folderPath = "C:/Users/Karthik/Documents/GitRepositories/Calculations/Imaging-Response-Function-Extractor/"; + +run = '0096'; + +folderPath = strcat(folderPath, run); + +cam = 5; + +angle = 0; +center = [1137, 2023]; +span = [500, 500]; +fraction = [0.1, 0.1]; + +pixel_size = 5.86e-6; +removeFringes = false; % Compute OD image, rotate and extract ROI for analysis % Get a list of all files in the folder with the desired file name pattern. @@ -36,8 +54,8 @@ for k = 1 : length(files) bkg_img = im2double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle)); dark_img = im2double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle)); - refimages(:,:,k) = subtract_offset(crop_image(bkg_img, center, span), fraction)'; - absimages(:,:,k) = subtract_offset(crop_image(calculate_OD(atm_img, bkg_img, dark_img), center, span), fraction)'; + refimages(:,:,k) = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)'; + absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img), center, span), fraction)'; end @@ -123,7 +141,7 @@ hold off; %% Helper Functions -function ret = get_offset_from_corner(img, x_fraction, y_fraction) +function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction) % image must be a 2D numerical array [dim1, dim2] = size(img); @@ -135,7 +153,7 @@ function ret = get_offset_from_corner(img, x_fraction, y_fraction) ret = mean([mean(s1(:)), mean(s2(:)), mean(s3(:)), mean(s4(:))]); end -function ret = subtract_offset(img, fraction) +function ret = subtractBackgroundOffset(img, fraction) % Remove the background from the image. % :param dataArray: The image % :type dataArray: xarray DataArray @@ -148,11 +166,11 @@ function ret = subtract_offset(img, fraction) x_fraction = fraction(1); y_fraction = fraction(2); - offset = get_offset_from_corner(img, x_fraction, y_fraction); + offset = getBkgOffsetFromCorners(img, x_fraction, y_fraction); ret = img - offset; end -function ret = crop_image(img, center, span) +function ret = cropODImage(img, center, span) % Crop the image according to the region of interest (ROI). % :param dataSet: The images % :type dataSet: xarray DataArray or DataSet @@ -171,7 +189,7 @@ function ret = crop_image(img, center, span) ret = img(y_start:y_end, x_start:x_end); end -function ret = calculate_OD(imageAtom, imageBackground, imageDark) +function ret = calculateODImage(imageAtom, imageBackground, imageDark) % Calculate the OD image for absorption imaging. % :param imageAtom: The image with atoms % :type imageAtom: numpy array @@ -264,4 +282,4 @@ function [optrefimages] = removefringesInImage(absimages, refimages, bgmask) % Compute optimised reference image optrefimages(:,:,j)=reshape(R*c,[ydim xdim]); end -end \ No newline at end of file +end diff --git a/Estimations/ModellingImagingAberrations.m b/Estimations/ModellingImagingAberrations.m new file mode 100644 index 0000000..b5bef97 --- /dev/null +++ b/Estimations/ModellingImagingAberrations.m @@ -0,0 +1,127 @@ +%% Zernike Polynomials + +resolution = 100; % Resolution + +plotZernike(2, 0, resolution) % Defocus (Z2^0) +%% +plotZernike(2, 2, resolution) % Astigmatism (Z2^2) +%% +plotZernike(3, 1, resolution) % Coma (Z3^1, x-direction) +%% +plotZernike(3, -1, resolution) % Coma (Z3^-1, y-direction) +%% +plotZernike(4, 0, resolution) % 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) + +plotAberratedPSF(C, pupil_radius, resolution); + +function Z = computeZernikePolynomials(n, m, r, theta) + % Zernike polynomial function for radial and angular coordinates (r, theta) + % Input: + % n - radial order + % m - azimuthal frequency + % r - radial coordinate (normalized to unit circle, 0 <= r <= 1) + % theta - angular coordinate (angle in radians) + % + % Output: + % Z - Zernike polynomial value at (r, theta) + + if n == 2 && m == 0 + % Defocus (Z2^0) + Z = 2 * r.^2 - 1; + elseif n == 2 && m == 2 + % Astigmatism (Z2^2) + Z = r.^2 .* cos(2 * theta); + elseif n == 3 && m == 1 + % Coma (Z3^1) + Z = (3 * r.^3 - 2 * r) .* cos(theta); + elseif n == 3 && m == -1 + % Coma (Z3^-1) + Z = (3 * r.^3 - 2 * r) .* sin(theta); + elseif n == 4 && m == 0 + % Spherical Aberration (Z4^0) + Z = 6 * r.^4 - 6 * r.^2 + 1; + else + % Default to zero if no known Zernike polynomial matches + Z = 0; + end +end + +function plotZernike(n, m, resolution) + % n: radial order + % m: azimuthal frequency + % resolution: number of points for plotting + + % Create a grid of (r, theta) coordinates + [theta, r] = meshgrid(linspace(0, 2*pi, resolution), linspace(0, 1, resolution)); + + % Calculate the Zernike polynomial for the given (n, m) + Z = computeZernikePolynomials(n, m, r, theta); + + % Convert polar to Cartesian coordinates for plotting + [X, Y] = pol2cart(theta, r); + + % Plot the Zernike polynomial using a surface plot + figure(1) + clf + set(gcf,'Position',[50 50 950 750]) + surf(X, Y, Z, 'EdgeColor', 'none'); + colormap jet; + colorbar; + title(sprintf('Zernike Polynomial Z_{%d}^{%d}', n, m), 'Interpreter', 'tex', 'FontSize', 16); + xlabel('X', 'FontSize', 16); + ylabel('Y', 'FontSize', 16); + zlabel('Zernike Value', 'FontSize', 16); + axis equal; + shading interp; +end + +function [theta, r, PSF] = modelPSF(C, pupil_radius, resolution) + % 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 + + % Create a grid of (r, theta) coordinates + [theta, r] = meshgrid(linspace(0, 2*pi, resolution), linspace(0, 1, resolution)); + + % Pupil function: 1 inside the pupil radius, 0 outside + P = double(r <= pupil_radius); % 2D mask + + % 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) + + % Fourier transform of the pupil function with aberrations + PSF = abs(fftshift(fft2(P .* exp(1i * W)))).^2; % Intensity distribution +end + +function plotAberratedPSF(C, pupil_radius, resolution) + % 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 + + % 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); + + figure(1) + clf + set(gcf,'Position',[50 50 950 750]) + surf(X, Y, PSF, 'EdgeColor', 'none'); + view(2); % 2D view + axis equal tight; + shading interp; + colorbar; + colormap jet; + title('PSF', 'FontSize', 16); + xlabel('X', 'FontSize', 16); + ylabel('Y', 'FontSize', 16); +end diff --git a/Imaging-Response-Function-Extractor/extractIRF.m b/Imaging-Response-Function-Extractor/extractIRF.m index 8845e61..2950cbe 100644 --- a/Imaging-Response-Function-Extractor/extractIRF.m +++ b/Imaging-Response-Function-Extractor/extractIRF.m @@ -1,26 +1,34 @@ %% Parameters -groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", "/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", "/images/Vertical_Axis_Camera/in_situ_absorption"]; +groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Camera/in_situ_absorption", "/images/ODT_2_Axis_Camera/in_situ_absorption", "/images/Horizontal_Axis_Camera/in_situ_absorption", "/images/Vertical_Axis_Camera/in_situ_absorption"]; -folderPath = "C:/Users/Karthik/Documents/GitRepositories/Calculations/Imaging-Response-Function-Extractor/0127/"; +folderPath = "C:/Users/Karthik/Documents/GitRepositories/Calculations/Imaging-Response-Function-Extractor/"; -cam = 5; +run = '0096'; -% angle = 90 + 51.5; -% center = [1700, 2300]; -angle = 90; -center = [2035 1250]; -span = [255, 255]; -fraction = [0.1, 0.1]; +folderPath = strcat(folderPath, run); -NA = 0.6; -pixel_size = 4.6e-6; -lambda = 421e-9; +cam = 5; -d = 0.61*lambda/NA; -k_cutoff = 2*NA/lambda/1e6; % in units of 1/µm) +angle = 0; +center = [1137, 2023]; +span = [255, 255]; +fraction = [0.1, 0.1]; -removeFringes = true; +NA = 0.6; +pixel_size = 4.6e-6; +lambda = 421e-9; + +% The diameter of the first dark concentric ring surrounding the central intensity peak of a point spread function (or Airy disk) +d = 1.22 * (lambda / NA); + +% Abbe limit +AbbeLimit = lambda / (2 * NA); + +% Maximum resolvable spatial frequency for the coherent case +k_cutoff = (NA/lambda) * 1e-6; % (in units of 1/µm) + +removeFringes = false; %% Compute OD image, rotate and extract ROI for analysis % Get a list of all files in the folder with the desired file name pattern. @@ -40,12 +48,12 @@ for k = 1 : length(files) bkg_img = im2double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle)); dark_img = im2double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle)); - refimages(:,:,k) = subtract_offset(crop_image(bkg_img, center, span), fraction); - absimages(:,:,k) = subtract_offset(crop_image(calculate_OD(atm_img, bkg_img, dark_img), center, span), fraction); + refimages(:,:,k) = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)'; + absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img), center, span), fraction)'; end -%% Fringe removal +% Fringe removal if removeFringes optrefimages = removefringesInImage(absimages, refimages); @@ -79,116 +87,109 @@ dy = pixel_size; xvals = (1:Nx)*dx*1e6; yvals = (1:Ny)*dy*1e6; -Nyq_k = 1/dx; % Nyquist -dk = 1/(Nx*dx); % Wavenumber increment -kx = -Nyq_k/2:dk:Nyq_k/2-dk; % wavenumber -kx = kx * dx; % wavenumber (in units of 1/dx) +dvx = 1 / (Nx * dx); % reciprocal space increment in the X direction (in units of 1/dx) +dvy = 1 / (Ny * dy); % reciprocal space increment in the Y direction (in units of 1/dy) -Nyq_k = 1/dy; % Nyquist -dk = 1/(Ny*dy); % Wavenumber increment -ky = -Nyq_k/2:dk:Nyq_k/2-dk; % wavenumber -ky = ky * dy; % wavenumber (in units of 1/dy) +% Create the frequency axes +vx = (-Nx/2:Nx/2-1) * dvx; % Frequency axis in X (in units of 1/dx) +vy = (-Ny/2:Ny/2-1) * dvy; % Frequency axis in Y (in units of 1/dy) + +% Create the Wavenumber axes +kx = 2*pi*vx; % Wavenumber axis in X +ky = 2*pi*vy; % Wavenumber axis in Y % Create Circular Mask -n = 2^8; % size of mask +n = 2^8; % size of mask mask = zeros(n); I = 1:n; -x = I-n/2; % mask x-coordinates -y = n/2-I; % mask y-coordinates -[X,Y] = meshgrid(x,y); % create 2-D mask grid -R = 32; % aperture radius -A = (X.^2 + Y.^2 <= R^2); % circular aperture of radius R -mask(A) = 1; % set mask elements inside aperture to 1 +x = I-n/2; % mask x-coordinates +y = n/2-I; % mask y-coordinates +[X,Y] = meshgrid(x,y); % create 2-D mask grid +R = 32; % aperture radius +A = (X.^2 + Y.^2 <= R^2); % circular aperture of radius R +mask(A) = 1; % set mask elements inside aperture to 1 % Calculate Power Spectrum and plot -figure('Position', [100, 100, 1200, 800]); +figure(1) clf +set(gcf,'Position',[100, 100, 1200, 800]) + +% Create tiled layout with 2 rows and 3 columns +t = tiledlayout(2, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); for k = 1 : length(od_imgs) mean_subtracted_od_imgs{k} = od_imgs{k} - mean_od_img; masked_img = mean_subtracted_od_imgs{k} .* mask; density_fft{k} = (1/numel(masked_img)) * abs(fftshift(fft2(masked_img))); density_noise_spectrum{k} = density_fft{k}.^2; - - % Subplot 1 - % subplot(2, 3, 1); - subplot('Position', [0.05, 0.55, 0.28, 0.4]) + + % Tile 1: Single-shot image + nexttile(1); imagesc(xvals, yvals, od_imgs{k}) - xlabel('µm', 'FontSize', 16) - ylabel('µm', 'FontSize', 16) + xlabel('\mum', 'Interpreter', 'tex', 'FontSize', 16) + ylabel('\mum', 'Interpreter', 'tex', 'FontSize', 16) axis equal tight; colorbar - colormap jet; % (flip(jet)) - % set(gca,'CLim',[0 100]); + colormap jet; set(gca,'YDir','normal') title('Single-shot image', 'FontSize', 16); - % Subplot 2 - % subplot(2, 3, 2); - subplot('Position', [0.36, 0.55, 0.28, 0.4]) + % Tile 2: Averaged density image + nexttile(2); imagesc(xvals, yvals, mean_od_img) - xlabel('µm', 'FontSize', 16) - ylabel('µm', 'FontSize', 16) + xlabel('\mum', 'Interpreter', 'tex', 'FontSize', 16) + ylabel('\mum', 'Interpreter', 'tex', 'FontSize', 16) axis equal tight; colorbar - colormap jet; % (flip(jet)) - % set(gca,'CLim',[0 100]); + colormap jet; set(gca,'YDir','normal') title('Averaged density image', 'FontSize', 16); - % Subplot 3 - % subplot(2, 3, 3); - subplot('Position', [0.67, 0.55, 0.28, 0.4]); + % Tile 3: Image noise = Single-shot - Average + nexttile(3); imagesc(xvals, yvals, mean_subtracted_od_imgs{k}) - xlabel('µm', 'FontSize', 16) - ylabel('µm', 'FontSize', 16) + xlabel('\mum', 'Interpreter', 'tex', 'FontSize', 16) + ylabel('\mum', 'Interpreter', 'tex', 'FontSize', 16) axis equal tight; colorbar - colormap jet; % (flip(jet)) - % set(gca,'CLim',[0 100]); + colormap jet; set(gca,'YDir','normal') title('Image noise = Single-shot - Average', 'FontSize', 16); - % Subplot 4 - % subplot(2, 3, 4); - subplot('Position', [0.05, 0.05, 0.28, 0.4]); + % Tile 4: Masked Noise + nexttile(4); imagesc(xvals, yvals, mean_subtracted_od_imgs{k} .* mask) - xlabel('µm', 'FontSize', 16) - ylabel('µm', 'FontSize', 16) + xlabel('\mum', 'Interpreter', 'tex', 'FontSize', 16) + ylabel('\mum', 'Interpreter', 'tex', 'FontSize', 16) axis equal tight; colorbar - colormap jet; % (flip(jet)) - % set(gca,'CLim',[0 100]); + colormap jet; set(gca,'YDir','normal') title('Masked Noise', 'FontSize', 16); - % Subplot 5 - % subplot(2, 3, 5); - subplot('Position', [0.36, 0.05, 0.28, 0.4]); - imagesc(kx, ky, abs(log2(density_fft{k}))) - xlabel('1/dx', 'FontSize', 16) - ylabel('1/dy', 'FontSize', 16) + % Tile 5: DFT + nexttile(5); + imagesc(kx*1E-6, ky*1E-6, abs(log2(density_fft{k}))) + xlabel('k_x (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16) + ylabel('k_y (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16) axis equal tight; colorbar - colormap jet; % (flip(jet)) - % set(gca,'CLim',[0 100]); + colormap jet; set(gca,'YDir','normal') title('DFT', 'FontSize', 16); - % Subplot 6 - % subplot(2, 3, 6); - subplot('Position', [0.67, 0.05, 0.28, 0.4]); - imagesc(kx, ky, abs(log2(density_noise_spectrum{k}))) - xlabel('1/dx', 'FontSize', 16) - ylabel('1/dy', 'FontSize', 16) + % Tile 6: Density Noise Spectrum = |DFT|^2 + nexttile(6); + imagesc(kx*1E-6, ky*1E-6, abs(log2(density_noise_spectrum{k}))) + xlabel('k_x (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16) + ylabel('k_y (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16) axis equal tight; colorbar - colormap jet; % (flip(jet)) - % set(gca,'CLim',[0 100]); + colormap jet; set(gca,'YDir','normal') - title('Density Noise Spectrum = |DFT|^2', 'FontSize', 16); - + title('Density Noise Spectrum', 'FontSize', 16); + drawnow; end @@ -198,25 +199,27 @@ end averagePowerSpectrum = mean(cat(3, density_noise_spectrum{:}), 3, 'double'); % Plot the average power spectrum. -figure('Position', [100, 100, 1200, 500]); +figure(2) clf +set(gcf,'Position',[100, 100, 1500, 700]) -subplot('Position', [0.05, 0.1, 0.4, 0.8]) % Adjusted position +% Create tiled layout with 2 rows and 3 columns +t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); + +nexttile(1); imagesc(abs(10*log10(averagePowerSpectrum))) axis equal tight; colorbar colormap(flip(jet)); -% set(gca,'CLim',[0 1e-7]); title('Average Density Noise Spectrum', 'FontSize', 16); grid on; + centers = ginput; -radius = 6; +radius = 3; % Plot where clicked. hVC = viscircles(centers, radius, 'Color', 'r', 'LineWidth', 2); xc = centers(:,1); -% xc = [78.2600, 108.3400, 128.8200, 150.5800, 181.3000]; yc = centers(:,2); -% yc = [131.3800, 155.7000, 128.8200, 101.3000, 126.2600]; [yDim, xDim] = size(averagePowerSpectrum); [xx,yy] = meshgrid(1:yDim,1:xDim); mask = false(xDim,yDim); @@ -225,78 +228,74 @@ for ii = 1:length(centers) end mask = not(mask); -x1 = 1; -y1 = 1; -x2 = 256; -y2 = 256; - % Ask user if the circle is acceptable. message = sprintf('Is this acceptable?'); button = questdlg(message, message, 'Accept', 'Reject and Quit', 'Accept'); if contains(button, 'Accept','IgnoreCase',true) image = mask.*averagePowerSpectrum; image(image==0) = NaN; - imagesc(kx, ky, mask.*abs(10*log10(averagePowerSpectrum))) + imagesc(kx*1E-6, ky*1E-6, mask.*abs(10*log10(averagePowerSpectrum))) hold on - line([kx(x1),kx(x2)], [ky(y1),ky(y2)], 'Color','white', 'LineStyle','--', 'LineWidth', 4); - % imagesc(kx, ky, 10*log10(averagePowerSpectrum)) - % imagesc(kx, ky, log2(averagePowerSpectrum)) - % imagesc(kx, ky, averagePowerSpectrum) - xlabel('1/dx', 'FontSize', 16) - ylabel('1/dy', 'FontSize', 16) + xlabel('k_x (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16) + ylabel('k_y (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16) axis equal tight; colorbar colormap(flip(jet)); - % set(gca,'CLim',[0 1e-7]); title('Average Density Noise Spectrum', 'FontSize', 16); grid on; elseif contains(button, 'Quit','IgnoreCase',true) delete(hVC); % Delete the circle from the overlay. image = averagePowerSpectrum; - imagesc(kx, ky, abs(10*log10(averagePowerSpectrum))) - % imagesc(kx, ky, 10*log10(averagePowerSpectrum)) - % imagesc(kx, ky, log2(averagePowerSpectrum)) - % imagesc(kx, ky, averagePowerSpectrum) - xlabel('1/dx', 'FontSize', 16) - ylabel('1/dy', 'FontSize', 16) + imagesc(kx*1E-6, ky*1E-6, abs(10*log10(averagePowerSpectrum))) + xlabel('k_x (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16) + ylabel('k_y (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16) axis equal tight; colorbar colormap(flip(jet)); - % set(gca,'CLim',[0 1e-7]); title('Average Density Noise Spectrum', 'FontSize', 16); grid on; end -subplot('Position', [0.55, 0.1, 0.4, 0.8]) % Adjusted position -% [r, Zr] = radial_profile(averagePowerSpectrum, 1); -% Zr = (Zr - min(Zr))./(max(Zr) - min(Zr)); -% plot(r, Zr, 'o-', 'MarkerSize', 4, 'MarkerFaceColor', 'none'); -% set(gca, 'XScale', 'log'); % Setting x-axis to log scale +% Fit +nexttile(2); +radialStep = 1; % in pixels +[R, Profile] = getRadialProfile(averagePowerSpectrum, radialStep); +kvec = (2 * pi) .* R(2:end) .* (sqrt(dvx^2 + dvy^2) * radialStep) * 1E-6; % in units of micrometers^-1 +NormalisedProfile = (Profile(2:end) - min(Profile(2:end)))./(max(Profile(2:end)) - min(Profile(2:end))); +kmax = k_cutoff; -[xi, yi, profile] = improfile(image, [x1,x2], [y1,y2]); -profile = (profile - min(profile))./(max(profile) - min(profile)); -ks = sqrt(kx.^2 + ky.^2); +% Define the objective function to minimize (the difference between model and data) +objectiveFunction = @(C, k) RadialImagingResponseFunction(C, k, kmax) - NormalisedProfile; -profile = profile(length(profile)/2:end); -ks = ks(length(ks)/2:end); +% Initial guess for the parameters [C1, C2, C3, C4, C5, C6] +initialGuess = [2E6, 1E-6, 1E-6, 1E-6, 1E-6, 0.8]; -n = 0.05; -[val,slice_idx]=min(abs(ks-n)); -ks = ks(1:slice_idx); -profile = profile(1:slice_idx); -plot(ks, profile, 'b*-'); -% plot(profile, 'b*-'); -grid on; -% xlim([min(ks) max(ks)]) -xlabel('k (1/µm)', 'FontSize', 16) +% Set upper and lower bounds for the parameters (optional) +lb = [-Inf, -Inf, -Inf, -Inf, -Inf, 0]; % Lower bounds +ub = [Inf, Inf, Inf, Inf, Inf, 1]; % Upper bounds + +% Perform the non-linear least squares fitting using lsqcurvefit +options = optimoptions('lsqcurvefit', 'Display', 'iter'); % Display iterations during fitting +[C_fit, resnorm] = lsqcurvefit(objectiveFunction, initialGuess, kvec, NormalisedProfile, lb, ub, options); + +% Plot the fitted result against the data +k_new = linspace(kvec(1), kvec(end), 1000); +fittedProfile = RadialImagingResponseFunction(C_fit, k_new, kmax); + +nexttile(2) +plot(kvec, NormalisedProfile, 'o-', 'MarkerSize', 4, 'MarkerFaceColor', 'none', 'DisplayName', 'Radial (average) profile'); +hold on +plot(k_new, fittedProfile, 'r-', 'DisplayName', 'Fitted Curve'); +set(gca, 'XScale', 'log'); % Setting x-axis to log scale +xlabel('k_\rho (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16) ylabel('Normalised amplitude', 'FontSize', 16) -title('Radial profile', 'FontSize', 16); +title('Modulation Transfer Function', 'FontSize', 16); +legend('FontSize', 16); grid on; - %% Helper Functions -function ret = get_offset_from_corner(img, x_fraction, y_fraction) +function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction) % image must be a 2D numerical array [dim1, dim2] = size(img); @@ -308,7 +307,7 @@ function ret = get_offset_from_corner(img, x_fraction, y_fraction) ret = mean([mean(s1(:)), mean(s2(:)), mean(s3(:)), mean(s4(:))]); end -function ret = subtract_offset(img, fraction) +function ret = subtractBackgroundOffset(img, fraction) % Remove the background from the image. % :param dataArray: The image % :type dataArray: xarray DataArray @@ -321,11 +320,11 @@ function ret = subtract_offset(img, fraction) x_fraction = fraction(1); y_fraction = fraction(2); - offset = get_offset_from_corner(img, x_fraction, y_fraction); + offset = getBkgOffsetFromCorners(img, x_fraction, y_fraction); ret = img - offset; end -function ret = crop_image(img, center, span) +function ret = cropODImage(img, center, span) % Crop the image according to the region of interest (ROI). % :param dataSet: The images % :type dataSet: xarray DataArray or DataSet @@ -344,7 +343,7 @@ function ret = crop_image(img, center, span) ret = img(y_start:y_end, x_start:x_end); end -function ret = calculate_OD(imageAtom, imageBackground, imageDark) +function ret = calculateODImage(imageAtom, imageBackground, imageDark) % Calculate the OD image for absorption imaging. % :param imageAtom: The image with atoms % :type imageAtom: numpy array @@ -368,50 +367,6 @@ function ret = calculate_OD(imageAtom, imageBackground, imageDark) end end -function [R, Zr] = radial_profile(data,radial_step) - x = (1:size(data,2))-size(data,2)/2; - y = (1:size(data,1))-size(data,1)/2; - % coordinate grid: - [X,Y] = meshgrid(x,y); - % creating circular layers - Z_integer = round(abs(X+1i*Y)/radial_step)+1; - % very fast MatLab calculations: - R = accumarray(Z_integer(:),abs(X(:)+1i*Y(:)),[],@mean); - Zr = accumarray(Z_integer(:),data(:),[],@mean); -end - -function [M] = ImagingResponseFunction(B) - x = -100:100; - y = x; - [X,Y] = meshgrid(x,y); - R = sqrt(X.^2+Y.^2); - PHI = atan2(X,Y)+pi; - %fit parameters - tau = B(1); - alpha = B(2); - S0 = B(3); - phi = B(4); - beta = B(5); - delta = B(6); - A = B(7); - C = B(8); - a = B(9); - U = heaviside(1-R/a).*exp(-R.^2/a^2/tau^2); - THETA = S0*(R/a).^4 + alpha*(R/a).^2.*cos(2*PHI-2*phi) + beta*(R/a).^2; - p = U.*exp(1i.*THETA); - M = A*abs((ifft2(real(exp(1i*delta).*fftshift(fft2(p)))))).^2 + C; -end - -function [RadialResponseFunc] = RadialImagingResponseFunction(C, k, kmax) - A = heaviside(1-k/kmax).*exp(-C(1)*k.^4); - W = C(2) + C(3)*k.^2 + C(4)*k.^4; - RadialResponseFunc = 0; - for n = -30:30 - RadialResponseFunc = RadialResponseFunc + besselj(n,C(5)*k.^2).^2 + besselj(n,C(5)*k.^2).*besselj(-n,C(5)*k.^2).*cos(2*W); - end - RadialResponseFunc = C(6)*1/2*A.*RadialResponseFunc; -end - function [optrefimages] = removefringesInImage(absimages, refimages, bgmask) % removefringesInImage - Fringe removal and noise reduction from absorption images. % Creates an optimal reference image for each absorption image in a set as @@ -481,4 +436,50 @@ function [optrefimages] = removefringesInImage(absimages, refimages, bgmask) % Compute optimised reference image optrefimages(:,:,j)=reshape(R*c,[ydim xdim]); end -end \ No newline at end of file +end + +function [M] = ImagingResponseFunction(B) + x = -100:100; + y = x; + [X,Y] = meshgrid(x,y); + R = sqrt(X.^2+Y.^2); + PHI = atan2(X,Y)+pi; + %fit parameters + tau = B(1); + alpha = B(2); + S0 = B(3); + phi = B(4); + beta = B(5); + delta = B(6); + A = B(7); + C = B(8); + a = B(9); + U = heaviside(1-R/a).*exp(-R.^2/a^2/tau^2); + THETA = S0*(R/a).^4 + alpha*(R/a).^2.*cos(2*PHI-2*phi) + beta*(R/a).^2; + p = U.*exp(1i.*THETA); + M = A*abs((ifft2(real(exp(1i*delta).*fftshift(fft2(p)))))).^2 + C; +end + +function [R, Zr] = getRadialProfile(data, radialStep) + x = (1:size(data,2))-size(data,2)/2; + y = (1:size(data,1))-size(data,1)/2; + % coordinate grid: + [X,Y] = meshgrid(x,y); + % creating circular layers + Z_integer = round(abs(X+1i*Y)/radialStep)+1; + % very fast MatLab calculations: + R = accumarray(Z_integer(:),abs(X(:)+1i*Y(:)),[],@mean); + Zr = accumarray(Z_integer(:),data(:),[],@mean); +end + +function [RadialResponseFunc] = RadialImagingResponseFunction(C, k, kmax) + A = heaviside(1 - k/kmax) .* exp(-C(1) * k.^4); + W = C(2) + C(3) * k.^2 + C(4) * k.^4; + RadialResponseFunc = 0; + for n = -30:30 + RadialResponseFunc = RadialResponseFunc + ... + besselj(n, C(5) * k.^2).^2 + ... + besselj(n, C(5) * k.^2) .* besselj(-n, C(5) * k.^2) .* cos(2 * W); + end + RadialResponseFunc = C(6) * 1/2 * A .* RadialResponseFunc; +end