Added direct correlation extraction and analysis of the OD images.

This commit is contained in:
Karthik 2025-09-15 11:38:44 +02:00
parent 20a0953b01
commit fd42630628
8 changed files with 544 additions and 4 deletions

View File

@ -0,0 +1,202 @@
function results = conductCorrelationAnalysis(od_imgs, scan_parameter_values, options)
%% conductCorrelationAnalysis
% Author: Karthik
% Date: 2025-09-14
% Version: 1.0
%
% Description:
% Computes 2D autocorrelation g²(Δx,Δy) on OD images.
% Extracts radial and angular distributions using Calculator functions.
% Optionally plots results and saves figures.
%
% Inputs:
% od_imgs - cell array of OD images
% scan_parameter_values - array of scan parameter values
% options - struct with fields:
% saveDirectory - base directory to save results
% skipSaveFigures - skip saving plots
% skipLivePlot - skip live plotting
% pixel_size - physical pixel size of camera sensor (m)
% magnification - imaging magnification
% maximumShift - maximum pixel shift for g²
% font - font name for plots
%
% Outputs:
% results - struct containing g² maps, radial and angular distributions
%
% Notes:
% Optional notes, references.
%% ===== Unpack struct arguments =====
pixel_size = options.pixel_size;
magnification = options.magnification;
skipLivePlot = options.skipLivePlot;
skipSaveFigures = options.skipSaveFigures;
saveDirectory = options.saveDirectory;
font = options.font;
radial_theta = options.Radial_Theta;
radial_window_size = options.Radial_WindowSize;
N_angular_bins = options.N_angular_bins;
r_min = options.Radial_Minimum;
r_max = options.Radial_Maximum;
if isfield(options, 'maximumShift') && ~isempty(options.maximumShift)
maximumShift = options.maximumShift;
else
maximumShift = 5; % [µm]
end
% --- Handle units ---
if ischar(options.scanParameterUnits) || isstring(options.scanParameterUnits)
unitList = {char(options.scanParameterUnits)};
else
unitList = options.scanParameterUnits;
end
%% ===== Initialization =====
N_shots = numel(od_imgs);
g2_matrices = cell(1, N_shots);
g2_radial = cell(1, N_shots);
g2_angular = cell(1, N_shots);
r_vals = cell(1, N_shots);
theta_vals = cell(1, N_shots);
dx = pixel_size / magnification;
shifts = -maximumShift:maximumShift;
if ~skipSaveFigures
saveFolder = fullfile(saveDirectory, 'Results', 'SavedFigures', 'AutocorrAnalysis');
if ~exist(saveFolder, 'dir')
mkdir(saveFolder);
end
end
%% ===== Loop over images =====
for k = 1:N_shots
IMG = od_imgs{k};
% Compute g² in physical units
[g2_matrix, dx_phys, dy_phys] = Calculator.compute2DAutocorrelation( ...
IMG, maximumShift, pixel_size, magnification);
g2_matrices{k} = g2_matrix;
% Extract radial profile (within angular window)
[r_vals{k}, g2_radial{k}] = Calculator.computeRadialCorrelation( ...
g2_matrix, dx_phys, dy_phys, radial_theta);
% Extract angular profile (within radial band)
[theta_vals{k}, g2_angular{k}] = Calculator.computeAngularCorrelation( ...
g2_matrix, dx_phys, dy_phys, r_min, r_max, N_angular_bins);
% Smooth radial profile
g2_radial_smoothed = movmean(g2_radial{k}, radial_window_size);
%% ===== Plotting =====
if ~skipLivePlot
figure(1); clf
set(gcf,'Position',[500 100 1000 800])
tiledlayout(2,2,'TileSpacing','compact','Padding','compact');
% OD image
ax1 = nexttile;
imagesc(IMG);
axis equal tight;
set(gca,'FontSize',14,'YDir','normal');
colormap(ax1, Colormaps.inferno());
hcb = colorbar;
ylabel(hcb,'Optical Density','Rotation',-90,'FontSize',14,'FontName',font);
xlabel('x [px]','FontSize',14,'FontName',font);
ylabel('y [px]','FontSize',14,'FontName',font);
title('OD Image','FontSize',16,'FontWeight','bold','FontName',font);
% Annotate scan parameter
if iscell(scan_parameter_values)
param_row = scan_parameter_values{k};
else
param_row = scan_parameter_values(k,:);
end
if numel(unitList) < numel(param_row)
unitList(end+1:numel(param_row)) = {''};
end
xPos = 0.975; yPos = 0.975; yStep = 0.075;
for j = 1:numel(param_row)
[unitSuffix, txtInterpreter] = getUnitInfo(unitList{j});
text(xPos, yPos-(j-1)*yStep, sprintf('%.2f%s',param_row(j),unitSuffix), ...
'Color','white','FontWeight','bold','FontSize',14, ...
'Interpreter',txtInterpreter,'Units','normalized', ...
'HorizontalAlignment','right','VerticalAlignment','top');
end
% g² map
ax2 = nexttile;
imagesc(shifts, shifts, g2_matrix);
axis equal tight;
set(gca,'FontSize',14,'YDir','normal');
colormap(ax2, Colormaps.coolwarm());
colorbar;
xlabel('\Deltax (\mum)','FontSize',14,'FontName',font);
ylabel('\Deltay (\mum)','FontSize',14,'FontName',font);
title('Autocorrelation g_2(\Deltax,\Deltay)','FontSize',16,'FontWeight','bold','FontName',font);
% Radial distribution
ax3 = nexttile;
plot(r_vals{k}, g2_radial_smoothed, 'LineWidth',2);
set(gca,'FontSize',14);
set(gca, 'FontSize', 14, 'XLim', [min(r_vals{k}), max(r_vals{k})], 'YLim', [0, 1]);
xlabel('r [\mum]','Interpreter','tex','FontSize',14,'FontName',font);
ylabel('g_2(r)','FontSize',14,'FontName',font);
title('Radial Distribution','FontSize',16,'FontWeight','bold','FontName',font);
grid on;
% Angular distribution
ax4 = nexttile;
plot(theta_vals{k}/pi, g2_angular{k}, 'LineWidth',2);
set(gca,'FontSize', 14, 'YLim', [0, 1]);
xlabel('\theta/\pi [rad]','Interpreter','tex','FontSize',14,'FontName',font);
ylabel('g_2(\theta)','FontSize',14,'FontName',font);
title('Angular Distribution','FontSize',16,'FontWeight','bold','FontName',font);
grid on;
ax = gca;
ax.MinorGridLineStyle = ':';
ax.MinorGridColor = [0.7 0.7 0.7];
ax.MinorGridAlpha = 0.5;
ax.XMinorGrid = 'on';
ax.YMinorGrid = 'on';
end
%% ===== Save figures =====
if ~skipSaveFigures && ~skipLivePlot
fileNamePNG = fullfile(saveFolder, sprintf('g2_analysis_img_%03d.png', k));
print(gcf, fileNamePNG, '-dpng','-r100');
elseif ~skipLivePlot
pause(0.5);
end
end
%% ===== Package results =====
results = struct();
results.g2_matrices = g2_matrices;
results.g2_radial = g2_radial;
results.g2_angular = g2_angular;
results.r_vals = r_vals;
results.theta_vals = theta_vals;
results.scan_parameter_values = unique(scan_parameter_values);
end
%% === Local helper function ===
function [unitSuffix, txtInterpreter] = getUnitInfo(u)
switch lower(u)
case {'degrees','deg','°'}
unitSuffix = '^\circ';
txtInterpreter = 'tex';
case {'gauss','g'}
unitSuffix = ' G';
txtInterpreter = 'none';
otherwise
unitSuffix = '';
txtInterpreter = 'none';
end
end

View File

@ -30,7 +30,7 @@ function results = conductSpectralAnalysis(od_imgs, scan_parameter_values, optio
% Angular_Sigma - Gaussian smoothing width for angular spectrum
% theta_min, theta_max - angular range for radial spectrum integration
% N_radial_bins - number of radial bins for S(k)
% Radial_WindowSize - window size for smoothing radial spectrum
% radial_window_size - window size for smoothing radial spectrum
% font - font name for plots
%
% Outputs:
@ -53,7 +53,7 @@ function results = conductSpectralAnalysis(od_imgs, scan_parameter_values, optio
theta_min = options.theta_min;
theta_max = options.theta_max;
N_radial_bins = options.N_radial_bins;
Radial_WindowSize = options.Radial_WindowSize;
radial_window_size = options.Radial_WindowSize;
skipNormalization = options.skipNormalization;
skipPreprocessing = options.skipPreprocessing;
skipMasking = options.skipMasking;
@ -144,7 +144,7 @@ function results = conductSpectralAnalysis(od_imgs, scan_parameter_values, optio
[k_rho_vals, S_k] = Calculator.computeRadialSpectralDistribution(fft_imgs{k}, kx, ky, theta_min, theta_max, N_radial_bins);
% Smooth radial spectrum
S_k_smoothed = movmean(S_k, Radial_WindowSize);
S_k_smoothed = movmean(S_k, radial_window_size);
% Store results
angular_spectral_distribution{k} = S_theta;

View File

@ -0,0 +1,87 @@
function results = extract2DAutocorrelation(od_imgs, scan_parameter_values, options)
%% extract2DAutocorrelation
% Author: Karthik
% Date: 2025-09-14
% Version: 1.0
%
% Description:
% Computes 2D autocorrelation g² for a set of real-space images.
% Returns all g2 maps grouped by unique scan parameter, mean, error
%
% Inputs:
% od_imgs - 1xN cell array of 2D image matrices (i × j)
% scan_parameter_values - array or cell array of scan parameters (similar format as in extractAutocorrelation)
% options - options struct containing maximum pixel shift to compute g2 along x and y
%
% Outputs (struct):
% results.g2_curves - cell array of 3D autocorrelation matrices [2*max_shift+1 × 2*max_shift+1 × N_reps] per unique scan parameter
% results.g2_mean - mean autocorrelation per unique scan parameter
% results.g2_error - SEM per unique scan parameter
% results.scan_parameter_values - unique scan parameter values
%
% Notes:
% Optional notes, references.
if isfield(options, 'maximumShift') && ~isempty(options.maximumShift)
max_shift = options.maximumShift;
else
max_shift = 5; % [µm]
end
pixel_size = options.pixel_size;
magnification = options.magnification;
% ===== Convert images to 3D array =====
N_shots = numel(od_imgs);
[M, N] = size(od_imgs{1});
img_stack = zeros(M, N, N_shots);
for k = 1:N_shots
img_stack(:, :, k) = od_imgs{k};
end
% ===== Determine unique scan parameter values =====
if isnumeric(scan_parameter_values) && isvector(scan_parameter_values)
[unique_scan_parameter_values, ~, idx] = unique(scan_parameter_values(:), 'stable');
elseif iscell(scan_parameter_values)
params = cell2mat(scan_parameter_values);
[unique_scan_parameter_values, ~, idx] = unique(params, 'rows', 'stable');
else
error('Unsupported format for scan_parameter_values.');
end
N_unique = size(unique_scan_parameter_values, 1);
% ===== Preallocate outputs =====
g2_curves = cell(1, N_unique); % each cell: [2*max_shift+1 × 2*max_shift+1 × N_reps]
% ===== Compute 2D autocorrelation for each unique scan parameter =====
for i = 1:N_unique
group_idx = find(idx == i); % indices of repetitions for this unique value
group_data = img_stack(:, :, group_idx);
N_reps = length(group_idx);
g2_stack = zeros(2*max_shift+1, 2*max_shift+1, N_reps);
for j = 1:N_reps
img = group_data(:, :, j);
[g2_matrix, ~, ~] = Calculator.compute2DAutocorrelation(img, max_shift, pixel_size, magnification);
g2_stack(:, :, j) = g2_matrix;
end
% Store all repetitions for this unique scan parameter
g2_curves{i} = g2_stack;
end
% ===== Compute mean and SEM per unique value =====
g2_mean = cellfun(@(G) mean(G, 3, 'omitnan'), g2_curves, 'UniformOutput', false);
g2_error = cellfun(@(G) std(G, 0, 3, 'omitnan') ./ sqrt(size(G,3)), g2_curves, 'UniformOutput', false);
% ===== Package results =====
results = struct();
results.g2_curves = g2_curves; % raw [2*max_shift+1 × 2*max_shift+1 × N_reps] per unique group
results.g2_mean = g2_mean; % mean per unique group
results.g2_error = g2_error; % SEM per unique group
results.scan_parameter_values = unique_scan_parameter_values;
results.max_shift = max_shift;
end

View File

@ -0,0 +1,49 @@
function [g2, dx_phys, dy_phys] = compute2DAutocorrelation(img, max_shift, pixel_size, magnification)
%% compute2DAutocorrelation
% Author: Karthik
% Date: 2025-09-14
% Version: 1.0
%
% Description:
% Computes 2D autocorrelation g2(Δx,Δy) with shifts in physical units [µm]
%
% Inputs:
% img : 2D image (OD)
% max_shift_um : maximum shift in micrometers
% pixel_size : camera pixel size in meters
% magnification: imaging magnification
%
% Outputs:
% g2 : normalized 2D autocorrelation
% dx_phys : x-axis shifts in µm
% dy_phys : y-axis shifts in µm
%
% Notes:
% Optional notes, references.
[M, N] = size(img);
img_mean_sub = img - mean(img(:));
I2_mean = mean(img_mean_sub(:).^2);
% Convert max_shift from µm to pixels
dx_max_px = round(max_shift / (pixel_size/magnification * 1e6));
g2 = zeros(2*dx_max_px+1, 2*dx_max_px+1);
for dx = -dx_max_px:dx_max_px
for dy = -dx_max_px:dx_max_px
% overlapping region in pixels
x1 = max(1,1+dx); x2 = min(M,M+dx);
y1 = max(1,1+dy); y2 = min(N,N+dy);
x1s = max(1,1-dx); x2s = min(M,M-dx);
y1s = max(1,1-dy); y2s = min(N,N-dy);
overlap = img_mean_sub(x1:x2, y1:y2) .* img_mean_sub(x1s:x2s, y1s:y2s);
g2(dx+dx_max_px+1, dy+dx_max_px+1) = mean(overlap(:)) / I2_mean;
end
end
% Return physical shift axes in µm
dx_phys = (-dx_max_px:dx_max_px) * (pixel_size/magnification * 1e6);
dy_phys = dx_phys;
end

View File

@ -0,0 +1,56 @@
function [theta_vals, g2_angular] = computeAngularCorrelation(g2_matrix, dx_phys, dy_phys, r_min, r_max, num_bins)
%% computeAngularCorrelation
% Author: Karthik
% Date: 2025-09-14
% Version: 1.0
%
% Description:
% Extracts angular profile of g2(Δx,Δy) along a radial band [r_min, r_max]
% from 0 to 180 degrees.
%
% Inputs:
% g2_matrix : 2D autocorrelation matrix
% dx_phys : x-axis shifts in µm
% dy_phys : y-axis shifts in µm
% r_min : minimum radial distance (µm)
% r_max : maximum radial distance (µm)
% num_bins : number of angular bins
%
% Outputs:
% theta_vals : angular positions [radians]
% g2_angular : angular profile of g2
%
% Notes:
% Optional notes, references.
[X, Y] = meshgrid(dx_phys, dy_phys);
R = sqrt(X.^2 + Y.^2);
Theta = atan2(Y, X); % [-pi, pi]
% Restrict to the radial band
radial_mask = (R >= r_min) & (R <= r_max);
% Angular bins from 0 to pi (0°180°)
theta_vals = linspace(0, pi, num_bins);
g2_angular = zeros(1, num_bins);
for i = 1:num_bins
angle_start = theta_vals(i) - (theta_vals(2)-theta_vals(1))/2;
angle_end = theta_vals(i) + (theta_vals(2)-theta_vals(1))/2;
% Handle wrap-around at 0/pi
if angle_start < 0
angle_mask = (Theta >= 0 & Theta < angle_end) | (Theta >= (pi+angle_start) & Theta <= pi);
elseif angle_end > pi
angle_mask = (Theta >= angle_start & Theta <= pi) | (Theta >= 0 & Theta < (angle_end - pi));
else
angle_mask = (Theta >= angle_start) & (Theta < angle_end);
end
% Combine with radial mask
bin_mask = radial_mask & angle_mask;
% Sum or average within this angular bin
g2_angular(i) = mean(g2_matrix(bin_mask), 'omitnan');
end
end

View File

@ -0,0 +1,44 @@
function [r_vals, g2_profile] = computeRadialCorrelation(g2_matrix, dx_phys, dy_phys, theta0)
%% computeRadialCorrelation
% Author: Karthik
% Date: 2025-09-14
% Version: 1.0
%
% Description:
% Extracts the g2 profile along a specific radial direction (line) at angle theta0
%
% Inputs:
% g2_matrix : 2D autocorrelation matrix
% dx_phys : x-axis shifts in µm
% dy_phys : y-axis shifts in µm
% theta0_deg : angle of the radial direction (degrees)
%
% Outputs:
% r_vals : radial distances along the line (µm)
% g2_profile : g2 values along that line
%
% Notes:
% Optional notes, references.
% Create meshgrid of physical shifts
[X, Y] = meshgrid(dx_phys, dy_phys);
% Compute radial distance along line
R = sqrt(X.^2 + Y.^2);
% Compute angle at each point
Theta = atan2(Y, X);
% Mask for points along the chosen line
% Use a small tolerance to pick points close to the line
tol = 0.5 * mean(diff(dx_phys)); % half pixel tolerance
line_mask = abs(Theta - theta0) < tol;
% Extract g2 values along the line
r_vals = R(line_mask);
g2_profile = g2_matrix(line_mask);
% Sort by radial distance
[r_vals, sort_idx] = sort(r_vals);
g2_profile = g2_profile(sort_idx);
end

View File

@ -89,5 +89,5 @@ end
[options.selectedPath, options.folderPath] = Helper.selectDataSourcePath(dataSources, options);
[od_imgs, scan_parameter_values, scan_reference_values, file_list] = Helper.collectODImages(options);
%%
Analyzer.runInteractiveODImageViewer(od_imgs, scan_parameter_values, file_list, options);

View File

@ -0,0 +1,102 @@
%% ===== BEC-Droplets-Stripes Settings =====
% Specify data location to run analysis on
dataSources = {
struct('sequence', 'TwoDGas', ...
'date', '2025/06/23', ...
'runs', [300]) % specify run numbers as a string in "" or just as a numeric value
};
options = struct();
% File paths
options.baseDataFolder = '//DyLabNAS/Data';
options.FullODImagesFolder = 'E:/Data - Experiment/FullODImages/202506';
options.measurementName = 'DropletsToStripes';
scriptFullPath = mfilename('fullpath');
options.saveDirectory = fileparts(scriptFullPath);
% Camera / imaging settings
options.cam = 4; % 1 - ODT_1_Axis_Camera; 2 - ODT_2_Axis_Camera; 3 - Horizontal_Axis_Camera;, 4 - Vertical_Axis_Camera;
options.angle = 0; % angle by which image will be rotated
options.center = [1410, 2030];
options.span = [200, 200];
options.fraction = [0.1, 0.1];
options.pixel_size = 5.86e-6; % in meters
options.magnification = 23.94;
options.ImagingMode = 'HighIntensity';
options.PulseDuration = 5e-6; % in s
% Fourier analysis settings
options.theta_min = deg2rad(0);
options.theta_max = deg2rad(180);
options.N_radial_bins = 500;
options.Radial_Sigma = 2;
options.Radial_WindowSize = 5; % odd number
options.k_min = 1.2771; % μm¹
options.k_max = 2.5541; % μm¹
options.N_angular_bins = 180;
options.Angular_Threshold = 75;
options.Angular_Sigma = 2;
options.Angular_WindowSize = 5;
options.zoom_size = 50;
%
options.maximumShift = 8;
options.Radial_Theta = deg2rad(45);
options.Radial_Minimum = 2;
options.Radial_Maximum = 6;
options.skipLivePlot = false;
% Flags
options.skipUnshuffling = false;
options.skipNormalization = false;
options.skipFringeRemoval = true;
options.skipPreprocessing = true;
options.skipMasking = true;
options.skipIntensityThresholding = true;
options.skipBinarization = true;
options.skipFullODImagesFolderUse = false;
options.skipSaveData = false;
options.skipSaveFigures = true;
options.skipSaveProcessedOD = true;
options.skipLivePlot = false;
options.showProgressBar = true;
% Extras
options.font = 'Bahnschrift';
switch options.measurementName
case 'BECToDroplets'
options.scan_parameter = 'rot_mag_field';
options.flipSortOrder = true;
options.scanParameterUnits = 'gauss';
options.titleString = 'BEC to Droplets';
case 'BECToStripes'
options.scan_parameter = 'rot_mag_field';
options.flipSortOrder = true;
options.scanParameterUnits = 'gauss';
options.titleString = 'BEC to Stripes';
case 'DropletsToStripes'
options.scan_parameter = 'ps_rot_mag_fin_pol_angle';
options.flipSortOrder = false;
options.scanParameterUnits = 'degrees';
options.titleString = 'Droplets to Stripes';
case 'StripesToDroplets'
options.scan_parameter = 'ps_rot_mag_fin_pol_angle';
options.flipSortOrder = false;
options.scanParameterUnits = 'degrees';
options.titleString = 'Stripes to Droplets';
end
%% ===== Collect Images and Launch Viewer =====
[options.selectedPath, options.folderPath] = Helper.selectDataSourcePath(dataSources, options);
[od_imgs, scan_parameter_values, scan_reference_values, file_list] = Helper.collectODImages(options);
%% Conduct correlation analysis
results = Analyzer.conductCorrelationAnalysis(od_imgs, scan_parameter_values, options);