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

Latest analysis script for IRF, new script to model imaging aberrations, modified plotting script for data from experiment.

Browse Source
This commit is contained in:
Karthik 2025-02-22 00:07:28 +01:00
parent 88c0d18066
commit 1ac6182f31
3 changed files with 337 additions and 191 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
Data-Analyzer/plotImages.m
View File

@ -1,12 +1,12 @@
%% Parameters %% 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; cam = 4;
@ -16,7 +16,25 @@ span = [50, 50];
fraction = [0.1, 0.1]; fraction = [0.1, 0.1];
pixel_size = 5.86e-6; 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 % Compute OD image, rotate and extract ROI for analysis
% Get a list of all files in the folder with the desired file name pattern. % 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)); bkg_img = im2double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle));
dark_img = im2double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle)); dark_img = im2double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle));
refimages(:,:,k) = subtract_offset(crop_image(bkg_img, center, span), fraction)'; refimages(:,:,k) = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
absimages(:,:,k) = subtract_offset(crop_image(calculate_OD(atm_img, bkg_img, dark_img), center, span), fraction)'; absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img), center, span), fraction)';
end end
@ -123,7 +141,7 @@ hold off;
%% Helper Functions %% 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 % image must be a 2D numerical array
[dim1, dim2] = size(img); [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(:))]); ret = mean([mean(s1(:)), mean(s2(:)), mean(s3(:)), mean(s4(:))]);
end end
function ret = subtract_offset(img, fraction) function ret = subtractBackgroundOffset(img, fraction)
% Remove the background from the image. % Remove the background from the image.
% :param dataArray: The image % :param dataArray: The image
% :type dataArray: xarray DataArray % :type dataArray: xarray DataArray
@ -148,11 +166,11 @@ function ret = subtract_offset(img, fraction)
x_fraction = fraction(1); x_fraction = fraction(1);
y_fraction = fraction(2); y_fraction = fraction(2);
offset = get_offset_from_corner(img, x_fraction, y_fraction); offset = getBkgOffsetFromCorners(img, x_fraction, y_fraction);
ret = img - offset; ret = img - offset;
end end
function ret = crop_image(img, center, span) function ret = cropODImage(img, center, span)
% Crop the image according to the region of interest (ROI). % Crop the image according to the region of interest (ROI).
% :param dataSet: The images % :param dataSet: The images
% :type dataSet: xarray DataArray or DataSet % :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); ret = img(y_start:y_end, x_start:x_end);
end end
function ret = calculate_OD(imageAtom, imageBackground, imageDark) function ret = calculateODImage(imageAtom, imageBackground, imageDark)
% Calculate the OD image for absorption imaging. % Calculate the OD image for absorption imaging.
% :param imageAtom: The image with atoms % :param imageAtom: The image with atoms
% :type imageAtom: numpy array % :type imageAtom: numpy array

127
Estimations/ModellingImagingAberrations.m Normal file
View File

@ -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

351
Imaging-Response-Function-Extractor/extractIRF.m
View File

@ -1,26 +1,34 @@
%% Parameters %% 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; folderPath = strcat(folderPath, run);
% center = [1700, 2300];
angle = 90;
center = [2035 1250];
span = [255, 255];
fraction = [0.1, 0.1];
NA = 0.6; cam = 5;
pixel_size = 4.6e-6;
lambda = 421e-9;
d = 0.61*lambda/NA; angle = 0;
k_cutoff = 2*NA/lambda/1e6; % in units of 1/µm) 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 %% Compute OD image, rotate and extract ROI for analysis
% Get a list of all files in the folder with the desired file name pattern. % 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)); bkg_img = im2double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle));
dark_img = im2double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle)); dark_img = im2double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle));
refimages(:,:,k) = subtract_offset(crop_image(bkg_img, center, span), fraction); refimages(:,:,k) = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
absimages(:,:,k) = subtract_offset(crop_image(calculate_OD(atm_img, bkg_img, dark_img), center, span), fraction); absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img), center, span), fraction)';
end end
%% Fringe removal % Fringe removal
if removeFringes if removeFringes
optrefimages = removefringesInImage(absimages, refimages); optrefimages = removefringesInImage(absimages, refimages);
@ -79,31 +87,36 @@ dy = pixel_size;
xvals = (1:Nx)*dx*1e6; xvals = (1:Nx)*dx*1e6;
yvals = (1:Ny)*dy*1e6; yvals = (1:Ny)*dy*1e6;
Nyq_k = 1/dx; % Nyquist dvx = 1 / (Nx * dx); % reciprocal space increment in the X direction (in units of 1/dx)
dk = 1/(Nx*dx); % Wavenumber increment dvy = 1 / (Ny * dy); % reciprocal space increment in the Y direction (in units of 1/dy)
kx = -Nyq_k/2:dk:Nyq_k/2-dk; % wavenumber
kx = kx * dx; % wavenumber (in units of 1/dx)
Nyq_k = 1/dy; % Nyquist % Create the frequency axes
dk = 1/(Ny*dy); % Wavenumber increment vx = (-Nx/2:Nx/2-1) * dvx; % Frequency axis in X (in units of 1/dx)
ky = -Nyq_k/2:dk:Nyq_k/2-dk; % wavenumber vy = (-Ny/2:Ny/2-1) * dvy; % Frequency axis in Y (in units of 1/dy)
ky = ky * dy; % wavenumber (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 % Create Circular Mask
n = 2^8; % size of mask n = 2^8; % size of mask
mask = zeros(n); mask = zeros(n);
I = 1:n; I = 1:n;
x = I-n/2; % mask x-coordinates x = I-n/2; % mask x-coordinates
y = n/2-I; % mask y-coordinates y = n/2-I; % mask y-coordinates
[X,Y] = meshgrid(x,y); % create 2-D mask grid [X,Y] = meshgrid(x,y); % create 2-D mask grid
R = 32; % aperture radius R = 32; % aperture radius
A = (X.^2 + Y.^2 <= R^2); % circular aperture of radius R A = (X.^2 + Y.^2 <= R^2); % circular aperture of radius R
mask(A) = 1; % set mask elements inside aperture to 1 mask(A) = 1; % set mask elements inside aperture to 1
% Calculate Power Spectrum and plot % Calculate Power Spectrum and plot
figure('Position', [100, 100, 1200, 800]); figure(1)
clf 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) for k = 1 : length(od_imgs)
mean_subtracted_od_imgs{k} = od_imgs{k} - mean_od_img; mean_subtracted_od_imgs{k} = od_imgs{k} - mean_od_img;
@ -111,83 +124,71 @@ for k = 1 : length(od_imgs)
density_fft{k} = (1/numel(masked_img)) * abs(fftshift(fft2(masked_img))); density_fft{k} = (1/numel(masked_img)) * abs(fftshift(fft2(masked_img)));
density_noise_spectrum{k} = density_fft{k}.^2; density_noise_spectrum{k} = density_fft{k}.^2;
% Subplot 1 % Tile 1: Single-shot image
% subplot(2, 3, 1); nexttile(1);
subplot('Position', [0.05, 0.55, 0.28, 0.4])
imagesc(xvals, yvals, od_imgs{k}) imagesc(xvals, yvals, od_imgs{k})
xlabel('µm', 'FontSize', 16) xlabel('\mum', 'Interpreter', 'tex', 'FontSize', 16)
ylabel('µm', 'FontSize', 16) ylabel('\mum', 'Interpreter', 'tex', 'FontSize', 16)
axis equal tight; axis equal tight;
colorbar colorbar
colormap jet; % (flip(jet)) colormap jet;
% set(gca,'CLim',[0 100]);
set(gca,'YDir','normal') set(gca,'YDir','normal')
title('Single-shot image', 'FontSize', 16); title('Single-shot image', 'FontSize', 16);
% Subplot 2 % Tile 2: Averaged density image
% subplot(2, 3, 2); nexttile(2);
subplot('Position', [0.36, 0.55, 0.28, 0.4])
imagesc(xvals, yvals, mean_od_img) imagesc(xvals, yvals, mean_od_img)
xlabel('µm', 'FontSize', 16) xlabel('\mum', 'Interpreter', 'tex', 'FontSize', 16)
ylabel('µm', 'FontSize', 16) ylabel('\mum', 'Interpreter', 'tex', 'FontSize', 16)
axis equal tight; axis equal tight;
colorbar colorbar
colormap jet; % (flip(jet)) colormap jet;
% set(gca,'CLim',[0 100]);
set(gca,'YDir','normal') set(gca,'YDir','normal')
title('Averaged density image', 'FontSize', 16); title('Averaged density image', 'FontSize', 16);
% Subplot 3 % Tile 3: Image noise = Single-shot - Average
% subplot(2, 3, 3); nexttile(3);
subplot('Position', [0.67, 0.55, 0.28, 0.4]);
imagesc(xvals, yvals, mean_subtracted_od_imgs{k}) imagesc(xvals, yvals, mean_subtracted_od_imgs{k})
xlabel('µm', 'FontSize', 16) xlabel('\mum', 'Interpreter', 'tex', 'FontSize', 16)
ylabel('µm', 'FontSize', 16) ylabel('\mum', 'Interpreter', 'tex', 'FontSize', 16)
axis equal tight; axis equal tight;
colorbar colorbar
colormap jet; % (flip(jet)) colormap jet;
% set(gca,'CLim',[0 100]);
set(gca,'YDir','normal') set(gca,'YDir','normal')
title('Image noise = Single-shot - Average', 'FontSize', 16); title('Image noise = Single-shot - Average', 'FontSize', 16);
% Subplot 4 % Tile 4: Masked Noise
% subplot(2, 3, 4); nexttile(4);
subplot('Position', [0.05, 0.05, 0.28, 0.4]);
imagesc(xvals, yvals, mean_subtracted_od_imgs{k} .* mask) imagesc(xvals, yvals, mean_subtracted_od_imgs{k} .* mask)
xlabel('µm', 'FontSize', 16) xlabel('\mum', 'Interpreter', 'tex', 'FontSize', 16)
ylabel('µm', 'FontSize', 16) ylabel('\mum', 'Interpreter', 'tex', 'FontSize', 16)
axis equal tight; axis equal tight;
colorbar colorbar
colormap jet; % (flip(jet)) colormap jet;
% set(gca,'CLim',[0 100]);
set(gca,'YDir','normal') set(gca,'YDir','normal')
title('Masked Noise', 'FontSize', 16); title('Masked Noise', 'FontSize', 16);
% Subplot 5 % Tile 5: DFT
% subplot(2, 3, 5); nexttile(5);
subplot('Position', [0.36, 0.05, 0.28, 0.4]); imagesc(kx*1E-6, ky*1E-6, abs(log2(density_fft{k})))
imagesc(kx, ky, abs(log2(density_fft{k}))) xlabel('k_x (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16)
xlabel('1/dx', 'FontSize', 16) ylabel('k_y (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16)
ylabel('1/dy', 'FontSize', 16)
axis equal tight; axis equal tight;
colorbar colorbar
colormap jet; % (flip(jet)) colormap jet;
% set(gca,'CLim',[0 100]);
set(gca,'YDir','normal') set(gca,'YDir','normal')
title('DFT', 'FontSize', 16); title('DFT', 'FontSize', 16);
% Subplot 6 % Tile 6: Density Noise Spectrum = |DFT|^2
% subplot(2, 3, 6); nexttile(6);
subplot('Position', [0.67, 0.05, 0.28, 0.4]); imagesc(kx*1E-6, ky*1E-6, abs(log2(density_noise_spectrum{k})))
imagesc(kx, ky, abs(log2(density_noise_spectrum{k}))) xlabel('k_x (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16)
xlabel('1/dx', 'FontSize', 16) ylabel('k_y (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16)
ylabel('1/dy', 'FontSize', 16)
axis equal tight; axis equal tight;
colorbar colorbar
colormap jet; % (flip(jet)) colormap jet;
% set(gca,'CLim',[0 100]);
set(gca,'YDir','normal') set(gca,'YDir','normal')
title('Density Noise Spectrum = |DFT|^2', 'FontSize', 16); title('Density Noise Spectrum', 'FontSize', 16);
drawnow; drawnow;
end end
@ -198,25 +199,27 @@ end
averagePowerSpectrum = mean(cat(3, density_noise_spectrum{:}), 3, 'double'); averagePowerSpectrum = mean(cat(3, density_noise_spectrum{:}), 3, 'double');
% Plot the average power spectrum. % Plot the average power spectrum.
figure('Position', [100, 100, 1200, 500]); figure(2)
clf 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))) imagesc(abs(10*log10(averagePowerSpectrum)))
axis equal tight; axis equal tight;
colorbar colorbar
colormap(flip(jet)); colormap(flip(jet));
% set(gca,'CLim',[0 1e-7]);
title('Average Density Noise Spectrum', 'FontSize', 16); title('Average Density Noise Spectrum', 'FontSize', 16);
grid on; grid on;
centers = ginput; centers = ginput;
radius = 6; radius = 3;
% Plot where clicked. % Plot where clicked.
hVC = viscircles(centers, radius, 'Color', 'r', 'LineWidth', 2); hVC = viscircles(centers, radius, 'Color', 'r', 'LineWidth', 2);
xc = centers(:,1); xc = centers(:,1);
% xc = [78.2600, 108.3400, 128.8200, 150.5800, 181.3000];
yc = centers(:,2); yc = centers(:,2);
% yc = [131.3800, 155.7000, 128.8200, 101.3000, 126.2600];
[yDim, xDim] = size(averagePowerSpectrum); [yDim, xDim] = size(averagePowerSpectrum);
[xx,yy] = meshgrid(1:yDim,1:xDim); [xx,yy] = meshgrid(1:yDim,1:xDim);
mask = false(xDim,yDim); mask = false(xDim,yDim);
@ -225,78 +228,74 @@ for ii = 1:length(centers)
end end
mask = not(mask); mask = not(mask);
x1 = 1;
y1 = 1;
x2 = 256;
y2 = 256;
% Ask user if the circle is acceptable. % Ask user if the circle is acceptable.
message = sprintf('Is this acceptable?'); message = sprintf('Is this acceptable?');
button = questdlg(message, message, 'Accept', 'Reject and Quit', 'Accept'); button = questdlg(message, message, 'Accept', 'Reject and Quit', 'Accept');
if contains(button, 'Accept','IgnoreCase',true) if contains(button, 'Accept','IgnoreCase',true)
image = mask.*averagePowerSpectrum; image = mask.*averagePowerSpectrum;
image(image==0) = NaN; 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 hold on
line([kx(x1),kx(x2)], [ky(y1),ky(y2)], 'Color','white', 'LineStyle','--', 'LineWidth', 4); xlabel('k_x (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16)
% imagesc(kx, ky, 10*log10(averagePowerSpectrum)) ylabel('k_y (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16)
% imagesc(kx, ky, log2(averagePowerSpectrum))
% imagesc(kx, ky, averagePowerSpectrum)
xlabel('1/dx', 'FontSize', 16)
ylabel('1/dy', 'FontSize', 16)
axis equal tight; axis equal tight;
colorbar colorbar
colormap(flip(jet)); colormap(flip(jet));
% set(gca,'CLim',[0 1e-7]);
title('Average Density Noise Spectrum', 'FontSize', 16); title('Average Density Noise Spectrum', 'FontSize', 16);
grid on; grid on;
elseif contains(button, 'Quit','IgnoreCase',true) elseif contains(button, 'Quit','IgnoreCase',true)
delete(hVC); % Delete the circle from the overlay. delete(hVC); % Delete the circle from the overlay.
image = averagePowerSpectrum; image = averagePowerSpectrum;
imagesc(kx, ky, abs(10*log10(averagePowerSpectrum))) imagesc(kx*1E-6, ky*1E-6, abs(10*log10(averagePowerSpectrum)))
% imagesc(kx, ky, 10*log10(averagePowerSpectrum)) xlabel('k_x (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16)
% imagesc(kx, ky, log2(averagePowerSpectrum)) ylabel('k_y (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16)
% imagesc(kx, ky, averagePowerSpectrum)
xlabel('1/dx', 'FontSize', 16)
ylabel('1/dy', 'FontSize', 16)
axis equal tight; axis equal tight;
colorbar colorbar
colormap(flip(jet)); colormap(flip(jet));
% set(gca,'CLim',[0 1e-7]);
title('Average Density Noise Spectrum', 'FontSize', 16); title('Average Density Noise Spectrum', 'FontSize', 16);
grid on; grid on;
end end
subplot('Position', [0.55, 0.1, 0.4, 0.8]) % Adjusted position % Fit
% [r, Zr] = radial_profile(averagePowerSpectrum, 1); nexttile(2);
% Zr = (Zr - min(Zr))./(max(Zr) - min(Zr)); radialStep = 1; % in pixels
% plot(r, Zr, 'o-', 'MarkerSize', 4, 'MarkerFaceColor', 'none'); [R, Profile] = getRadialProfile(averagePowerSpectrum, radialStep);
% set(gca, 'XScale', 'log'); % Setting x-axis to log scale 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]); % Define the objective function to minimize (the difference between model and data)
profile = (profile - min(profile))./(max(profile) - min(profile)); objectiveFunction = @(C, k) RadialImagingResponseFunction(C, k, kmax) - NormalisedProfile;
ks = sqrt(kx.^2 + ky.^2);
profile = profile(length(profile)/2:end); % Initial guess for the parameters [C1, C2, C3, C4, C5, C6]
ks = ks(length(ks)/2:end); initialGuess = [2E6, 1E-6, 1E-6, 1E-6, 1E-6, 0.8];
n = 0.05; % Set upper and lower bounds for the parameters (optional)
[val,slice_idx]=min(abs(ks-n)); lb = [-Inf, -Inf, -Inf, -Inf, -Inf, 0]; % Lower bounds
ks = ks(1:slice_idx); ub = [Inf, Inf, Inf, Inf, Inf, 1]; % Upper bounds
profile = profile(1:slice_idx);
plot(ks, profile, 'b*-'); % Perform the non-linear least squares fitting using lsqcurvefit
% plot(profile, 'b*-'); options = optimoptions('lsqcurvefit', 'Display', 'iter'); % Display iterations during fitting
grid on; [C_fit, resnorm] = lsqcurvefit(objectiveFunction, initialGuess, kvec, NormalisedProfile, lb, ub, options);
% xlim([min(ks) max(ks)])
xlabel('k (1/µm)', 'FontSize', 16) % 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) ylabel('Normalised amplitude', 'FontSize', 16)
title('Radial profile', 'FontSize', 16); title('Modulation Transfer Function', 'FontSize', 16);
legend('FontSize', 16);
grid on; grid on;
%% Helper Functions %% 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 % image must be a 2D numerical array
[dim1, dim2] = size(img); [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(:))]); ret = mean([mean(s1(:)), mean(s2(:)), mean(s3(:)), mean(s4(:))]);
end end
function ret = subtract_offset(img, fraction) function ret = subtractBackgroundOffset(img, fraction)
% Remove the background from the image. % Remove the background from the image.
% :param dataArray: The image % :param dataArray: The image
% :type dataArray: xarray DataArray % :type dataArray: xarray DataArray
@ -321,11 +320,11 @@ function ret = subtract_offset(img, fraction)
x_fraction = fraction(1); x_fraction = fraction(1);
y_fraction = fraction(2); y_fraction = fraction(2);
offset = get_offset_from_corner(img, x_fraction, y_fraction); offset = getBkgOffsetFromCorners(img, x_fraction, y_fraction);
ret = img - offset; ret = img - offset;
end end
function ret = crop_image(img, center, span) function ret = cropODImage(img, center, span)
% Crop the image according to the region of interest (ROI). % Crop the image according to the region of interest (ROI).
% :param dataSet: The images % :param dataSet: The images
% :type dataSet: xarray DataArray or DataSet % :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); ret = img(y_start:y_end, x_start:x_end);
end end
function ret = calculate_OD(imageAtom, imageBackground, imageDark) function ret = calculateODImage(imageAtom, imageBackground, imageDark)
% Calculate the OD image for absorption imaging. % Calculate the OD image for absorption imaging.
% :param imageAtom: The image with atoms % :param imageAtom: The image with atoms
% :type imageAtom: numpy array % :type imageAtom: numpy array
@ -368,50 +367,6 @@ function ret = calculate_OD(imageAtom, imageBackground, imageDark)
end end
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) function [optrefimages] = removefringesInImage(absimages, refimages, bgmask)
% removefringesInImage - Fringe removal and noise reduction from absorption images. % removefringesInImage - Fringe removal and noise reduction from absorption images.
% Creates an optimal reference image for each absorption image in a set as % Creates an optimal reference image for each absorption image in a set as
@ -482,3 +437,49 @@ function [optrefimages] = removefringesInImage(absimages, refimages, bgmask)
optrefimages(:,:,j)=reshape(R*c,[ydim xdim]); optrefimages(:,:,j)=reshape(R*c,[ydim xdim]);
end end
end 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