Calculations/Data-Analyzer/extractCustomCorrelation.m

646 lines
24 KiB
Matlab
Raw Blame History

This file contains invisible Unicode characters

This file contains invisible Unicode characters that are indistinguishable to humans but may be processed differently by a computer. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

%% ===== Settings =====
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 = "//DyLabNAS/Data/TwoDGas/2025/06/23/";
run = '0300';
folderPath = strcat(folderPath, run);
cam = 5;
angle = 0;
center = [1410, 2030];
span = [200, 200];
fraction = [0.1, 0.1];
pixel_size = 5.86e-6; % in meters
magnification = 23.94;
removeFringes = false;
ImagingMode = 'HighIntensity';
PulseDuration = 5e-6; % in s
% Fourier analysis settings
% Radial Spectral Distribution
theta_min = deg2rad(0);
theta_max = deg2rad(180);
N_radial_bins = 500;
Radial_Sigma = 2;
Radial_WindowSize = 5; % Choose an odd number for a centered moving average
% Angular Spectral Distribution
r_min = 10;
r_max = 20;
N_angular_bins = 180;
Angular_Threshold = 75;
Angular_Sigma = 2;
Angular_WindowSize = 5;
zoom_size = 50; % Zoomed-in region around center
% Plotting and saving
scan_parameter = 'ps_rot_mag_fin_pol_angle';
% scan_parameter = 'rot_mag_field';
scan_parameter_text = 'Angle = ';
% scan_parameter_text = 'BField = ';
savefileName = 'DropletsToStripes';
font = 'Bahnschrift';
if strcmp(savefileName, 'DropletsToStripes')
scan_groups = 0:5:45;
titleString = 'Droplets to Stripes';
elseif strcmp(savefileName, 'StripesToDroplets')
scan_groups = 45:-5:0;
titleString = 'Stripes to Droplets';
end
% Flags
skipUnshuffling = true;
skipPreprocessing = true;
skipMasking = true;
skipIntensityThresholding = true;
skipBinarization = true;
skipMovieRender = true;
skipSaveFigures = true;
%% ===== Load and compute OD image, rotate and extract ROI for analysis =====
% Get a list of all files in the folder with the desired file name pattern.
filePattern = fullfile(folderPath, '*.h5');
files = dir(filePattern);
refimages = zeros(span(1) + 1, span(2) + 1, length(files));
absimages = zeros(span(1) + 1, span(2) + 1, length(files));
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
fprintf(1, 'Now reading %s\n', fullFileName);
atm_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/atoms")), angle));
bkg_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/background")), angle));
dark_img = double(imrotate(h5read(fullFileName, append(groupList(cam), "/dark")), angle));
refimages(:,:,k) = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)';
absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)';
end
%% ===== Fringe removal =====
if removeFringes
optrefimages = removefringesInImage(absimages, refimages);
absimages_fringe_removed = absimages(:, :, :) - optrefimages(:, :, :);
nimgs = size(absimages_fringe_removed,3);
od_imgs = cell(1, nimgs);
for i = 1:nimgs
od_imgs{i} = absimages_fringe_removed(:, :, i);
end
else
nimgs = size(absimages(:, :, :),3);
od_imgs = cell(1, nimgs);
for i = 1:nimgs
od_imgs{i} = absimages(:, :, i);
end
end
%% ===== Get rotation angles =====
scan_parameter_values = zeros(1, length(files));
% Get information about the '/globals' group
for k = 1 : length(files)
baseFileName = files(k).name;
fullFileName = fullfile(files(k).folder, baseFileName);
info = h5info(fullFileName, '/globals');
for i = 1:length(info.Attributes)
if strcmp(info.Attributes(i).Name, scan_parameter)
if strcmp(scan_parameter, 'ps_rot_mag_fin_pol_angle')
scan_parameter_values(k) = 180 - info.Attributes(i).Value;
else
scan_parameter_values(k) = info.Attributes(i).Value;
end
end
end
end
%% ===== Correlation of a single (highest) peak with a possible peak between 50-70 degrees from experiment data =====
fft_imgs = cell(1, nimgs);
spectral_distribution = cell(1, nimgs);
theta_values = cell(1, nimgs);
N_shots = length(od_imgs);
% Compute FFT for all images
for k = 1:N_shots
IMG = od_imgs{k};
[IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization);
[Ny, Nx] = size(IMG);
dx = pixel_size / magnification;
dy = dx; % assuming square pixels
x = ((1:Nx) - ceil(Nx/2)) * dx * 1E6;
y = ((1:Ny) - ceil(Ny/2)) * dy * 1E6;
dvx = 1 / (Nx * dx);
dvy = 1 / (Ny * dy);
vx = (-floor(Nx/2):ceil(Nx/2)-1) * dvx;
vy = (-floor(Ny/2):ceil(Ny/2)-1) * dvy;
kx_full = 2 * pi * vx * 1E-6;
ky_full = 2 * pi * vy * 1E-6;
mid_x = floor(Nx/2);
mid_y = floor(Ny/2);
fft_imgs{k} = IMGFFT(mid_y-zoom_size:mid_y+zoom_size, mid_x-zoom_size:mid_x+zoom_size);
kx = kx_full(mid_x - zoom_size : mid_x + zoom_size);
ky = ky_full(mid_y - zoom_size : mid_y + zoom_size);
[theta_vals, S_theta] = computeAngularSpectralDistribution(fft_imgs{k}, r_min, r_max, N_angular_bins, Angular_Threshold, Angular_Sigma, []);
spectral_distribution{k} = S_theta;
theta_values{k} = theta_vals;
end
% Convert spectral distribution to matrix (N_shots x N_angular_bins)
delta_nkr_all = zeros(N_shots, N_angular_bins);
for k = 1:N_shots
delta_nkr_all(k, :) = spectral_distribution{k};
end
% Group by scan parameter values (e.g., alpha, angle, etc.)
[unique_scan_parameter_values, ~, idx] = unique(scan_parameter_values);
N_params = length(unique_scan_parameter_values);
% Define angular range and conversion
angle_range = 180;
angle_per_bin = angle_range / N_angular_bins;
max_peak_angle = 180;
max_peak_bin = round(max_peak_angle / angle_per_bin);
% Parameters for search
window_size = 10;
angle_threshold = 100;
% Initialize containers for final results
mean_max_g2_values = zeros(1, N_params);
mean_max_g2_angle_values = zeros(1, N_params);
var_max_g2_values = zeros(1, N_params);
var_max_g2_angle_values = zeros(1, N_params);
% Also store raw data per group
g2_all_per_group = cell(1, N_params);
angle_all_per_group = cell(1, N_params);
for i = 1:N_params
group_idx = find(idx == i);
group_data = delta_nkr_all(group_idx, :);
N_reps = size(group_data, 1);
g2_values = zeros(1, N_reps);
angle_at_max_g2 = zeros(1, N_reps);
for j = 1:N_reps
profile = group_data(j, :);
% Restrict search to 0–60° for highest peak
restricted_profile = profile(1:max_peak_bin);
[~, peak_idx_rel] = max(restricted_profile);
peak_idx = peak_idx_rel;
peak_angle = (peak_idx - 1) * angle_per_bin;
if peak_angle < angle_threshold
offsets = round(50 / angle_per_bin) : round(70 / angle_per_bin);
else
offsets = -round(70 / angle_per_bin) : -round(50 / angle_per_bin);
end
ref_window = mod((peak_idx - window_size):(peak_idx + window_size) - 1, N_angular_bins) + 1;
ref = profile(ref_window);
correlations = zeros(size(offsets));
angles = zeros(size(offsets));
for k = 1:length(offsets)
shifted_idx = mod(peak_idx + offsets(k) - 1, N_angular_bins) + 1;
sec_window = mod((shifted_idx - window_size):(shifted_idx + window_size) - 1, N_angular_bins) + 1;
sec = profile(sec_window);
num = mean(ref .* sec);
denom = mean(ref.^2);
g2 = num / denom;
correlations(k) = g2;
angles(k) = mod((peak_idx - 1 + offsets(k)) * angle_per_bin, angle_range);
end
[max_corr, max_idx] = max(correlations);
g2_values(j) = max_corr;
angle_at_max_g2(j) = angles(max_idx);
end
% Store raw values
g2_all_per_group{i} = g2_values;
angle_all_per_group{i} = angle_at_max_g2;
% Final stats
mean_max_g2_values(i) = mean(g2_values, 'omitnan');
var_max_g2_values(i) = var(g2_values, 0, 'omitnan');
mean_max_g2_angle_values(i)= mean(angle_at_max_g2, 'omitnan');
var_max_g2_angle_values(i) = var(angle_at_max_g2, 0, 'omitnan');
end
%% ── Mean ± Std vs. scan parameter ──────────────────────────────────────
% Compute standard error instead of standard deviation
std_error_g2_values = zeros(1, N_params);
for i = 1:N_params
n_i = numel(g2_all_per_group{i}); % Number of repetitions for this param
std_error_g2_values(i) = sqrt(var_max_g2_values(i) / n_i);
end
% Plot mean ± SEM
figure(1);
set(gcf,'Position',[100 100 950 750])
set(gca, 'FontSize', 14); % For tick labels only
errorbar(unique_scan_parameter_values, ... % x-axis
mean_max_g2_values, ... % y-axis (mean)
std_error_g2_values, ... % ± SEM
'--o', 'LineWidth', 1.8, 'MarkerSize', 6 );
set(gca, 'FontSize', 14, 'YLim', [0, 1]);
hXLabel = xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex');
hYLabel = ylabel('$\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', 'Interpreter', 'latex');
hTitle = title(titleString, 'Interpreter', 'tex');
% set([hXLabel, hYLabel], 'FontName', font);
set([hXLabel, hYLabel], 'FontSize', 14);
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold');
grid on;
% Define folder for saving images
saveFolder = [savefileName '_SavedFigures'];
if ~exist(saveFolder, 'dir')
mkdir(saveFolder);
end
save([saveFolder savefileName '.mat'], 'unique_scan_parameter_values', 'mean_max_g2_values', 'std_error_g2_values');
%{
%% Plot histograms within 0-180 degrees only
figure(1);
hold on;
bin_edges = 0:10:180;
h1 = histogram(ref_peak_angles, 'BinEdges', bin_edges, ...
'FaceColor', [0.3 0.7 0.9], 'EdgeColor', 'none', 'FaceAlpha', 0.6);
h2 = histogram(angle_at_max_g2, 'BinEdges', bin_edges, ...
'FaceColor', [0.9 0.4 0.4], 'EdgeColor', 'none', 'FaceAlpha', 0.6);
h1.Normalization = 'probability';
h2.Normalization = 'probability';
xlabel('Angle (degrees)', 'FontSize', 12);
ylabel('Probability', 'FontSize', 12);
legend({'Reference Peak Angle', 'Angle at Max gâ‚‚'}, 'FontSize', 12);
title('Comparison of Reference Peak and Max gâ‚‚ Angles', 'FontSize', 14);
grid on;
xlim([0 180]);
hold off;
% Assume ref_peak_angles and angle_at_max_g2 are row or column vectors of angles in [0,180]
% Define fine angle grid for KDE evaluation
angle_grid = linspace(0, 180, 1000);
% KDE for reference peak angles
[f_ref, xi_ref] = ksdensity(ref_peak_angles, angle_grid, 'Bandwidth', 5);
% KDE for max g2 angles
[f_g2, xi_g2] = ksdensity(angle_at_max_g2, angle_grid, 'Bandwidth', 5);
% Plot KDEs
figure(2);
plot(xi_ref, f_ref, 'LineWidth', 2, 'DisplayName', 'Reference Peak Angles');
hold on;
plot(xi_g2, f_g2, 'LineWidth', 2, 'DisplayName', 'Max g_2 Angles');
xlabel('Angle (degrees)');
ylabel('Probability Density');
title('KDE of Angle Distributions');
legend;
grid on;
% Find modes (angle at max KDE value)
[~, mode_idx_ref] = max(f_ref);
mode_ref = xi_ref(mode_idx_ref);
[~, mode_idx_g2] = max(f_g2);
mode_g2 = xi_g2(mode_idx_g2);
% Calculate difference in mode
mode_diff = abs(mode_ref - mode_g2);
fprintf('Mode difference between distributions: %.2f degrees\n', mode_diff);
% Add vertical dashed lines at mode positions
yl = ylim; % get y-axis limits for text positioning
% Reference peak mode line and label
xline(mode_ref, 'k--', 'LineWidth', 1.5, 'DisplayName', sprintf('Ref Mode: %.1f°', mode_ref));
text(mode_ref, yl(2)*0.9, sprintf('%.1f°', mode_ref), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 12, 'Color', 'k');
% Max g2 mode line and label
xline(mode_g2, 'r--', 'LineWidth', 1.5, 'DisplayName', sprintf('g_2 Mode: %.1f°', mode_g2));
text(mode_g2, yl(2)*0.75, sprintf('%.1f°', mode_g2), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 12, 'Color', 'r');
%}
%% Helper Functions
function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization)
% computeFourierSpectrum - Computes the 2D Fourier power spectrum
% of binarized and enhanced lattice image features, with optional central mask.
%
% Inputs:
% I - Grayscale or RGB image matrix
%
% Output:
% F_mag - 2D Fourier power spectrum (shifted)
if ~skipPreprocessing
% Preprocessing: Denoise
filtered = imgaussfilt(I, 10);
IMGPR = I - filtered; % adjust sigma as needed
else
IMGPR = I;
end
if ~skipMasking
[rows, cols] = size(IMGPR);
[X, Y] = meshgrid(1:cols, 1:rows);
% Elliptical mask parameters
cx = cols / 2;
cy = rows / 2;
% Shifted coordinates
x = X - cx;
y = Y - cy;
% Ellipse semi-axes
rx = 0.4 * cols;
ry = 0.2 * rows;
% Rotation angle in degrees -> radians
theta_deg = 30; % Adjust as needed
theta = deg2rad(theta_deg);
% Rotated ellipse equation
cos_t = cos(theta);
sin_t = sin(theta);
x_rot = (x * cos_t + y * sin_t);
y_rot = (-x * sin_t + y * cos_t);
ellipseMask = (x_rot.^2) / rx^2 + (y_rot.^2) / ry^2 <= 1;
% Apply cutout mask
IMGPR = IMGPR .* ellipseMask;
end
if ~skipIntensityThresholding
% Apply global intensity threshold mask
intensity_thresh = 0.20;
intensity_mask = IMGPR > intensity_thresh;
IMGPR = IMGPR .* intensity_mask;
end
if ~skipBinarization
% Adaptive binarization and cleanup
IMGPR = imbinarize(IMGPR, 'adaptive', 'Sensitivity', 0.0);
IMGPR = imdilate(IMGPR, strel('disk', 2));
IMGPR = imerode(IMGPR, strel('disk', 1));
IMGPR = imfill(IMGPR, 'holes');
F = fft2(double(IMGPR)); % Compute 2D Fourier Transform
IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center
else
F = fft2(double(IMGPR)); % Compute 2D Fourier Transform
IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center
end
end
function [theta_vals, S_theta] = computeAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma, windowSize)
% Apply threshold to isolate strong peaks
IMGFFT(IMGFFT < threshold) = 0;
% Prepare polar coordinates
[ny, nx] = size(IMGFFT);
[X, Y] = meshgrid(1:nx, 1:ny);
cx = ceil(nx/2);
cy = ceil(ny/2);
R = sqrt((X - cx).^2 + (Y - cy).^2);
Theta = atan2(Y - cy, X - cx); % range [-pi, pi]
% Choose radial band
radial_mask = (R >= r_min) & (R <= r_max);
% Initialize angular structure factor
S_theta = zeros(1, num_bins);
theta_vals = linspace(0, pi, num_bins);
% Loop through angle bins
for i = 1:num_bins
angle_start = (i-1) * pi / num_bins;
angle_end = i * pi / num_bins;
angle_mask = (Theta >= angle_start & Theta < angle_end);
bin_mask = radial_mask & angle_mask;
fft_angle = IMGFFT .* bin_mask;
S_theta(i) = sum(sum(abs(fft_angle).^2));
end
% Smooth using either Gaussian or moving average
if exist('sigma', 'var') && ~isempty(sigma)
% Gaussian convolution
half_width = ceil(3 * sigma);
x = -half_width:half_width;
gauss_kernel = exp(-x.^2 / (2 * sigma^2));
gauss_kernel = gauss_kernel / sum(gauss_kernel);
% Circular convolution
S_theta = conv([S_theta(end-half_width+1:end), S_theta, S_theta(1:half_width)], ...
gauss_kernel, 'same');
S_theta = S_theta(half_width+1:end-half_width);
elseif exist('windowSize', 'var') && ~isempty(windowSize)
% Moving average via convolution (circular)
pad = floor(windowSize / 2);
kernel = ones(1, windowSize) / windowSize;
S_theta = conv([S_theta(end-pad+1:end), S_theta, S_theta(1:pad)], kernel, 'same');
S_theta = S_theta(pad+1:end-pad);
end
end
function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction)
% image must be a 2D numerical array
[dim1, dim2] = size(img);
s1 = img(1:round(dim1 * y_fraction), 1:round(dim2 * x_fraction));
s2 = img(1:round(dim1 * y_fraction), round(dim2 - dim2 * x_fraction):dim2);
s3 = img(round(dim1 - dim1 * y_fraction):dim1, 1:round(dim2 * x_fraction));
s4 = img(round(dim1 - dim1 * y_fraction):dim1, round(dim2 - dim2 * x_fraction):dim2);
ret = mean([mean(s1(:)), mean(s2(:)), mean(s3(:)), mean(s4(:))]);
end
function ret = subtractBackgroundOffset(img, fraction)
% Remove the background from the image.
% :param dataArray: The image
% :type dataArray: xarray DataArray
% :param x_fraction: The fraction of the pixels used in x axis
% :type x_fraction: float
% :param y_fraction: The fraction of the pixels used in y axis
% :type y_fraction: float
% :return: The image after removing background
% :rtype: xarray DataArray
x_fraction = fraction(1);
y_fraction = fraction(2);
offset = getBkgOffsetFromCorners(img, x_fraction, y_fraction);
ret = img - offset;
end
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
% :param center: The center of region of interest (ROI)
% :type center: tuple
% :param span: The span of region of interest (ROI)
% :type span: tuple
% :return: The cropped images
% :rtype: xarray DataArray or DataSet
x_start = floor(center(1) - span(1) / 2);
x_end = floor(center(1) + span(1) / 2);
y_start = floor(center(2) - span(2) / 2);
y_end = floor(center(2) + span(2) / 2);
ret = img(y_start:y_end, x_start:x_end);
end
function imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%CALCULATEODIMAGE Calculates the optical density (OD) image for absorption imaging.
%
% imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime)
%
% Inputs:
% imageAtom - Image with atoms
% imageBackground - Image without atoms
% imageDark - Image without light
% mode - 'LowIntensity' (default) or 'HighIntensity'
% exposureTime - Required only for 'HighIntensity' [in seconds]
%
% Output:
% imageOD - Computed OD image
%
arguments
imageAtom (:,:) {mustBeNumeric}
imageBackground (:,:) {mustBeNumeric}
imageDark (:,:) {mustBeNumeric}
mode char {mustBeMember(mode, {'LowIntensity', 'HighIntensity'})} = 'LowIntensity'
exposureTime double = NaN
end
% Compute numerator and denominator
numerator = imageBackground - imageDark;
denominator = imageAtom - imageDark;
% Avoid division by zero
numerator(numerator == 0) = 1;
denominator(denominator == 0) = 1;
% Calculate OD based on mode
switch mode
case 'LowIntensity'
imageOD = -log(abs(denominator ./ numerator));
case 'HighIntensity'
if isnan(exposureTime)
error('Exposure time must be provided for HighIntensity mode.');
end
imageOD = abs(denominator ./ numerator);
imageOD = -log(imageOD) + (numerator - denominator) ./ (7000 * (exposureTime / 5e-6));
end
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
% a linear combination of reference images, with coefficients chosen to
% minimize the least-squares residuals between each absorption image and
% the optimal reference image. The coefficients are obtained by solving a
% linear set of equations using matrix inverse by LU decomposition.
%
% Application of the algorithm is described in C. F. Ockeloen et al, Improved
% detection of small atom numbers through image processing, arXiv:1007.2136 (2010).
%
% Syntax:
% [optrefimages] = removefringesInImage(absimages,refimages,bgmask);
%
% Required inputs:
% absimages - Absorption image data,
% typically 16 bit grayscale images
% refimages - Raw reference image data
% absimages and refimages are both cell arrays containing
% 2D array data. The number of refimages can differ from the
% number of absimages.
%
% Optional inputs:
% bgmask - Array specifying background region used,
% 1=background, 0=data. Defaults to all ones.
% Outputs:
% optrefimages - Cell array of optimal reference images,
% equal in size to absimages.
%
% Dependencies: none
%
% Authors: Shannon Whitlock, Caspar Ockeloen
% Reference: C. F. Ockeloen, A. F. Tauschinsky, R. J. C. Spreeuw, and
% S. Whitlock, Improved detection of small atom numbers through
% image processing, arXiv:1007.2136
% Email:
% May 2009; Last revision: 11 August 2010
% Process inputs
% Set variables, and flatten absorption and reference images
nimgs = size(absimages,3);
nimgsR = size(refimages,3);
xdim = size(absimages(:,:,1),2);
ydim = size(absimages(:,:,1),1);
R = single(reshape(refimages,xdim*ydim,nimgsR));
A = single(reshape(absimages,xdim*ydim,nimgs));
optrefimages=zeros(size(absimages)); % preallocate
if not(exist('bgmask','var')); bgmask=ones(ydim,xdim); end
k = find(bgmask(:)==1); % Index k specifying background region
% Ensure there are no duplicate reference images
% R=unique(R','rows')'; % comment this line if you run out of memory
% Decompose B = R*R' using singular value or LU decomposition
[L,U,p] = lu(R(k,:)'*R(k,:),'vector'); % LU decomposition
for j=1:nimgs
b=R(k,:)'*A(k,j);
% Obtain coefficients c which minimise least-square residuals
lower.LT = true; upper.UT = true;
c = linsolve(U,linsolve(L,b(p,:),lower),upper);
% Compute optimised reference image
optrefimages(:,:,j)=reshape(R*c,[ydim xdim]);
end
end