diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzeFolder.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzeFolder.m index ed483b3..e0ac599 100644 --- a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzeFolder.m +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzeFolder.m @@ -96,8 +96,8 @@ function results = analyzeFolder(options) (isempty(bkg_img) && isa(bkg_img, 'double')) || ... (isempty(dark_img) && isa(dark_img, 'double')) - refimages = nan(size(refimages)); % fill with NaNs - absimages = nan(size(absimages)); + refimages(:,:,k) = nan(size(refimages(:,:,k))); % fill with NaNs + absimages(:,:,k) = nan(size(absimages(:,:,k))); else 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)'; diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzewithPCA.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzewithPCA.m index b741159..ead36e6 100644 --- a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzewithPCA.m +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/analyzewithPCA.m @@ -1,210 +1,505 @@ -%% STEP 1: Extract Images +%% Extract Images clear; close all; clc; -baseDir = 'D:/Results - Numerics/Data_Full3D/PhaseTransition/DTS/Point2'; -JobNumber = 0; -runFolder = sprintf('Run_%03d', JobNumber); -savefileName = 'DropletsToStripes'; % Output file name -datafileName = './DropletsToStripes.mat'; -reverseOrder = false; % Set this to true to reverse the theta ordering +%% ===== D-S 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"]; -TitleString = 'Droplets To Stripes'; +folderPath = "//DyLabNAS/Data/TwoDGas/2025/06/23/"; -%% -baseDir = 'D:/Results - Numerics/Data_Full3D/PhaseTransition/STD/Point2'; -JobNumber = 0; -runFolder = sprintf('Run_%03d', JobNumber); -savefileName = 'StripesToDroplets'; % Output file name -datafileName = './StripesToDroplets.mat'; -reverseOrder = true; % Set this to true to reverse the theta ordering +run = '0300'; -TitleString = 'Stripes To Droplets'; +folderPath = strcat(folderPath, run); -%% -folderList = dir(baseDir); -isValid = [folderList.isdir] & ~ismember({folderList.name}, {'.', '..'}); -folderNames = {folderList(isValid).name}; -nimgs = numel(folderNames); +cam = 5; -% Extract theta values from folder names -PolarAngleVals = zeros(1, nimgs); +angle = 0; +center = [1410, 2030]; +span = [200, 200]; +fraction = [0.1, 0.1]; -for k = 1:nimgs - tokens = regexp(folderNames{k}, 'theta_(\d{3})', 'tokens'); - if isempty(tokens) - warning('No theta found in folder name: %s', folderNames{k}); - PolarAngleVals(k) = NaN; - else - PolarAngleVals(k) = str2double(tokens{1}{1}); - end +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'; + +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 -% Choose sort direction -sortDirection = 'ascend'; -if reverseOrder - sortDirection = 'descend'; -end +% Flags +skipNormalization = true; +skipUnshuffling = true; +skipPreprocessing = true; +skipMasking = true; +skipIntensityThresholding = true; +skipBinarization = true; +skipMovieRender = true; +skipSaveFigures = true; +skipSaveOD = true; -% Sort folderNames based on polar angle -[~, sortIdx] = sort(PolarAngleVals, sortDirection); -folderNames = folderNames(sortIdx); -PolarAngleVals = PolarAngleVals(sortIdx); % Optional: if you still want sorted list -imgs = cell(1, nimgs); -alphas = zeros(1, nimgs); +%% ===== 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. -for k = 1:nimgs - folderName = folderNames{k}; - SaveDirectory = fullfile(baseDir, folderName, runFolder); +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)); - % Extract alpha (theta) again from folder name - tokens = regexp(folderName, 'theta_(\d{3})', 'tokens'); - alpha_val = str2double(tokens{1}{1}); - alphas(k) = alpha_val; - - matPath = fullfile(SaveDirectory, 'psi_gs.mat'); - if ~isfile(matPath) - warning('Missing psi_gs.mat in %s', SaveDirectory); - continue; - end - - try - Data = load(matPath, 'psi', 'Params', 'Transf', 'Observ'); - catch ME - warning('Failed to load %s: %s', matPath, ME.message); - continue; - end - - Params = Data.Params; - Transf = Data.Transf; - Observ = Data.Observ; - - psi = Data.psi; - if isgpuarray(psi) - psi = gather(psi); - end - if isgpuarray(Observ.residual) - Observ.residual = gather(Observ.residual); - end +for k = 1 : length(files) + baseFileName = files(k).name; + fullFileName = fullfile(files(k).folder, baseFileName); - % Axes and projection - z = Transf.z * Params.l0 * 1e6; - dz = z(2)-z(1); + fprintf(1, 'Now reading %s\n', fullFileName); - if k == 1 - x = Transf.x * Params.l0 * 1e6; - y = Transf.y * Params.l0 * 1e6; - % Calculate frequency increment (frequency axes) - Nx = length(x); % grid size along X - Ny = length(y); % grid size along Y - dx = mean(diff(x)); % real space increment in the X direction (in micrometers) - dy = mean(diff(y)); % real space increment in the Y direction (in micrometers) - dvx = 1 / (Nx * dx); % reciprocal space increment in the X direction (in micrometers^-1) - dvy = 1 / (Ny * dy); % reciprocal space increment in the Y direction (in micrometers^-1) + 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)); + + if (isempty(atm_img) && isa(atm_img, 'double')) || ... + (isempty(bkg_img) && isa(bkg_img, 'double')) || ... + (isempty(dark_img) && isa(dark_img, 'double')) - % Create the frequency axes - vx = (-Nx/2:Nx/2-1) * dvx; % Frequency axis in X (micrometers^-1) - vy = (-Ny/2:Ny/2-1) * dvy; % Frequency axis in Y (micrometers^-1) - - % Create the Wavenumber axes - kx_full = 2*pi*vx; % Wavenumber axis in X - ky_full = 2*pi*vy; % Wavenumber axis in Y + refimages(:,:,k) = nan(size(refimages(:,:,k))); % fill with NaNs + absimages(:,:,k) = nan(size(absimages(:,:,k))); + else + 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 - - n = abs(psi).^2; - nxy = squeeze(trapz(n * dz, 3)); - - imgs{k} = nxy; end -%% STEP 2: Visualize images -figure(1); -set(gcf,'Position',[100 100 950 750]) -font = 'Bahnschrift'; -t = tiledlayout(4, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid - -for i = 1:nimgs - ax1 = nexttile; - imagesc(x, y, imgs{i}') - % Define normalized positions (relative to axis limits) - x_offset = 0.025; % 5% offset from the edges - y_offset = 0.025; % 5% offset from the edges - axis square; - hcb = colorbar; - colormap(ax1, Colormaps.plasma()); - set(gca, 'FontSize', 14); % For tick labels only - hL = ylabel(hcb, 'Optical Density'); - set(hL,'Rotation',-90); - set(gca,'YDir','normal') - % set(gca, 'YTick', linspace(y_min, y_max, 5)); % Define y ticks - % set(gca, 'YTickLabel', flip(linspace(y_min, y_max, 5))); % Flip only the labels - hXLabel = xlabel('x (µm)', 'Interpreter', 'tex'); - hYLabel = ylabel('y (µm)', 'Interpreter', 'tex'); - hTitle = title('Density', 'Interpreter', 'tex'); - set([hXLabel, hYLabel, hL], 'FontName', font) - set([hXLabel, hYLabel, hL], 'FontSize', 14) - set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title +% ===== 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 -sgtitle('Ground States'); +% ===== Get rotation angles ===== +scan_parameter_values = zeros(1, length(files)); -%% STEP 3: Reshape images into a 2D data matrix -% Each image will be a row vector, so size(X) = [numImages × numPixels] +% 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 -% First, convert cell array to numeric 3D array -allImgs = cat(3, imgs{:}); % size: [Nx × Ny × numImages] +% ===== Unshuffle if necessary to do so ===== -% Store image size for reshaping later -imgSize = size(allImgs(:,:,1)); % [Nx, Ny] +if ~skipUnshuffling + n_values = length(scan_groups); + n_total = length(scan_parameter_values); + + % Infer number of repetitions + n_reps = n_total / n_values; + + % Preallocate ordered arrays + ordered_scan_values = zeros(1, n_total); + ordered_od_imgs = cell(1, n_total); + + counter = 1; + + for rep = 1:n_reps + for val = scan_groups + % Find the next unused match for this val + idx = find(scan_parameter_values == val, 1, 'first'); + + % Assign and remove from list to avoid duplicates + ordered_scan_values(counter) = scan_parameter_values(idx); + ordered_od_imgs{counter} = od_imgs{idx}; + + % Mark as used by removing + scan_parameter_values(idx) = NaN; % NaN is safe since original values are 0:5:45 + od_imgs{idx} = []; % empty cell so it won't be matched again + + counter = counter + 1; + end + end + + % Now assign back + scan_parameter_values = ordered_scan_values; + od_imgs = ordered_od_imgs; +end -% Reshape into 2D matrix for PCA -X = reshape(allImgs, [], nimgs)'; % [numImages × (Nx*Ny)] +%% Carry out PCA +numPCs = 5; -%% STEP 4: Perform PCA -[coeff, score, latent, tsquared, explained] = pca(X); +% Stack all 600 images into one data matrix [nImages x nPixels] +allImgs3D = cat(3, od_imgs{:}); +[Nx, Ny] = size(allImgs3D(:,:,1)); +Xall = reshape(allImgs3D, [], numel(od_imgs))'; % [600 x (Nx*Ny)] -% coeff : principal component directions (eigenvectors) in pixel space -% score : projection of each image onto the principal components -% latent : eigenvalues (variance captured by each PC) -% explained : percentage variance explained by each PC +% Global PCA +[coeff, score, ~, ~, explained] = pca(Xall); -%% Visualize cumulative variance -figure(2); -set(gcf,'Position',[100 100 950 750]) -cumVar = cumsum(explained); -plot(cumVar, 'o-', 'LineWidth', 1.5); -yline(95, '--r', '95% threshold'); -xlabel('Number of Principal Components'); -ylabel('Cumulative Variance (%)'); -title('Minimum number of PCA components required for a basis'); +%% Visualize PC1 +% Extract the first principal component vector (eigenimage) +pc1_vector = coeff(:,1); + +% Reshape back to original image dimensions +pc1_image = reshape(pc1_vector, Nx, Ny); + +% Plot the PC1 image +figure(1); clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]); +imagesc(pc1_image); +axis image off; +colormap(Colormaps.coolwarm()); % or use 'jet', 'parula', etc. +colorbar; +title(sprintf('First Principal Component (PC1) Image - Explains %.2f%% Variance', explained(1))); + +%% Distribution scatter plot +colors = lines(numGroups); + +figure(2); clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]); hold on; +for g = 1:numGroups + idx = scan_parameter_values == scan_groups(g); + scatter(repmat(scan_groups(g), sum(idx),1), score(idx,1), 36, colors(g,:), 'filled'); +end +xlabel('Control Parameter'); +ylabel('PC1 Score'); +title('Evolution of PC1 Scores'); grid on; -%% Visualize first few principal components -numPCsToShow = 7; -figure(3); -set(gcf,'Position',[100 100 950 750]) -t = tiledlayout(3, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid -for i = 1:numPCsToShow - pcImage = reshape(coeff(:,i), imgSize); % reshape back to image form - nexttile; - imagesc(pcImage) - axis image off; colormap(Colormaps.coolwarm()); - title(sprintf('PC %d (%.1f%%)', i, explained(i))); -end -sgtitle('First Principal Components (Eigenimages)'); +%% Distribution Histogram plot +numGroups = length(scan_groups); +colors = lines(numGroups); + +% Define number of bins globally +numBins = 20; + +% Define common bin edges based on global PC1 score range +minScore = min(score(:,1)); +maxScore = max(score(:,1)); +binEdges = linspace(minScore, maxScore, numBins+1); % +1 because edges are one more than bins +binWidth = binEdges(2) - binEdges(1); % for scaling KDE + +figure(3); +clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]); +tiledlayout(ceil(numGroups/2), 2, 'TileSpacing', 'compact', 'Padding', 'compact'); + +for g = 1:numGroups + groupVal = scan_groups(g); + idx = scan_parameter_values == groupVal; + groupPC1 = score(idx,1); + + nexttile; + + % Plot histogram + histogram(groupPC1, 'Normalization', 'probability', ... + 'FaceColor', colors(g,:), 'EdgeColor', 'none', ... + 'BinEdges', binEdges); + hold on; + + % Compute KDE + [f, xi] = ksdensity(groupPC1, 'NumPoints', 1000); + + % Scale KDE to histogram probability scale + f_scaled = f * binWidth; + + % Overlay KDE curve + plot(xi, f_scaled, 'k', 'LineWidth', 1.5); + + % Vertical line at median + med = median(groupPC1); + yl = ylim; + plot([med med], yl, 'k--', 'LineWidth', 1); + + xlabel('PC1 Score'); + ylabel('Probability'); + title(sprintf('Control Parameter = %d', groupVal)); + grid on; + hold off; +end + +sgtitle('PC1 Score Distributions'); + +%% Box plot for PC1 scores by group +groupLabels = cell(size(score,1),1); +for g = 1:numGroups + idx = scan_parameter_values == scan_groups(g); + groupLabels(idx) = {sprintf('%d', scan_groups(g))}; +end -%% Visualize evolution of PCA scores figure(4); -set(gcf,'Position',[100 100 950 750]) -plot(alphas, score(:,1), '-o', 'LineWidth', 1.5, 'DisplayName', 'PC1 Score'); -hold on; -plot(alphas, score(:,2), '-o', 'LineWidth', 1.5, 'DisplayName', 'PC2 Score'); -% plot(alphas, score(:,3), '-o', 'LineWidth', 1.5, 'DisplayName', 'PC3 Score'); -% plot(alphas, score(:,4), '-o', 'LineWidth', 1.5, 'DisplayName', 'PC4 Score'); -% plot(alphas, score(:,5), '-o', 'LineWidth', 1.5, 'DisplayName', 'PC5 Score'); -% plot(alphas, score(:,6), '-o', 'LineWidth', 1.5, 'DisplayName', 'PC6 Score'); -% plot(alphas, score(:,7), '-o', 'LineWidth', 1.5, 'DisplayName', 'PC7 Score'); -legend('Location', 'northeast'); -xlabel('Polar Angle (°)'); ylabel('PC Score'); -title('Evolution of PC Scores'); -grid on; \ No newline at end of file +clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]); +boxplot(score(:,1), groupLabels); +xlabel('Control Parameter'); +ylabel('PC1 Score'); +title('Evolution of PC1 Scores'); +grid on; + +%% Mean and SEM plot for PC1 scores +numGroups = length(scan_groups); +meanPC1Scores = zeros(numGroups,1); +semPC1Scores = zeros(numGroups,1); + +for g = 1:numGroups + groupVal = scan_groups(g); + idx = scan_parameter_values == groupVal; + groupPC1 = score(idx,1); % PC1 scores for this group + + meanPC1Scores(g) = mean(groupPC1); + semPC1Scores(g) = std(groupPC1)/sqrt(sum(idx)); % Standard error of mean +end + +% Plot mean ± SEM with error bars +figure(5); +clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]); +errorbar(scan_groups, meanPC1Scores, semPC1Scores, 'o-', ... + 'LineWidth', 1.5, 'MarkerSize', 8, 'MarkerFaceColor', 'b'); +xlabel('Control Parameter'); +ylabel('Mean PC1 Score ± SEM'); +title('Evolution of PC1 Scores'); +grid on; + +%% Plot Binder Cumulant +maxOrder = 4; % We only need up to order 4 here +numGroups = length(scan_groups); +kappa4 = NaN(1, numGroups); + +for g = 1:numGroups + groupVal = scan_groups(g); + idx = scan_parameter_values == groupVal; + groupPC1 = score(idx, 1); + + cumulants = computeCumulants(groupPC1, maxOrder); + kappa4(g) = cumulants(4); % 4th-order cumulant +end + +% Plot +figure(6); +clf; set(gcf, 'Color', 'w', 'Position', [100 100 950 750]); +plot(scan_groups, kappa4 * 1E-5, '-o', 'LineWidth', 1.5, 'MarkerFaceColor', 'b'); +ylim([-12 12]) +xlabel('Control Parameter'); +ylabel('\kappa_4 (\times 10^{5})'); +grid on; +title('Evolution of Binder Cumulant of PC1 Score'); + +%% --- ANOVA test --- +p = anova1(score(:,1), groupLabels, 'off'); +fprintf('ANOVA p-value for PC1 score differences between groups: %.4e\n', p); + +%% Helper Functions + +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 + diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/conductSpectralAnalysis.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/conductSpectralAnalysis.m index 11922ca..9c0b052 100644 --- a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/conductSpectralAnalysis.m +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/conductSpectralAnalysis.m @@ -3,9 +3,9 @@ groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1 "/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/07/22/"; +folderPath = "//DyLabNAS/Data/TwoDGas/2025/06/23/"; -run = '0021'; +run = '0300'; folderPath = strcat(folderPath, run); @@ -50,16 +50,16 @@ savefileName = 'DropletsToStripes'; font = 'Bahnschrift'; if strcmp(savefileName, 'DropletsToStripes') - scan_groups = 0:1:40; + scan_groups = 0:5:45; titleString = 'Droplets to Stripes'; elseif strcmp(savefileName, 'StripesToDroplets') - scan_groups = 40:-1:0; + scan_groups = 45:-5:0; titleString = 'Stripes to Droplets'; end % Flags skipNormalization = true; -skipUnshuffling = false; +skipUnshuffling = true; skipPreprocessing = true; skipMasking = true; skipIntensityThresholding = true; diff --git a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractQuantities.m b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractQuantities.m index dccb69f..98168aa 100644 --- a/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractQuantities.m +++ b/Data-Analyzer/StructuralPhaseTransition/SpectralAnalysisRoutines/extractQuantities.m @@ -10,29 +10,29 @@ runs = { ["0007", "0008", "0009", "0010", "0011"] }; -options.scan_parameter = 'rot_mag_fin_pol_angle'; -options.scan_groups = 0:10:50; -options.cam = 5; +options.scan_parameter = 'rot_mag_fin_pol_angle'; +options.scan_groups = 0:10:50; +options.cam = 5; % Image cropping and alignment -options.angle = 0; -options.center = [1285, 2100]; -options.span = [200, 200]; -options.fraction = [0.1, 0.1]; +options.angle = 0; +options.center = [1285, 2100]; +options.span = [200, 200]; +options.fraction = [0.1, 0.1]; % Imaging and calibration parameters -options.pixel_size = 5.86e-6; % in meters -options.magnification = 23.94; -options.removeFringes = false; -options.ImagingMode = 'HighIntensity'; -options.PulseDuration = 5e-6; +options.pixel_size = 5.86e-6; % in meters +options.magnification = 23.94; +options.removeFringes = false; +options.ImagingMode = 'HighIntensity'; +options.PulseDuration = 5e-6; % Fourier analysis: Radial -options.theta_min = deg2rad(0); -options.theta_max = deg2rad(180); -options.N_radial_bins = 500; -options.Radial_Sigma = 2; -options.Radial_WindowSize = 5; % Must be odd +options.theta_min = deg2rad(0); +options.theta_max = deg2rad(180); +options.N_radial_bins = 500; +options.Radial_Sigma = 2; +options.Radial_WindowSize = 5; % Must be odd % Fourier analysis: Angular options.r_min = 10; @@ -45,14 +45,14 @@ options.Angular_Sigma = 2; options.Angular_WindowSize = 5; % Optional visualization / zooming -options.zoom_size = 50; +options.zoom_size = 50; % Optional flags or settings struct -options.skipUnshuffling = false; -options.skipPreprocessing = true; -options.skipMasking = true; +options.skipUnshuffling = false; +options.skipPreprocessing = true; +options.skipMasking = true; options.skipIntensityThresholding = true; -options.skipBinarization = true; +options.skipBinarization = true; % === Loop through folders and collect results === @@ -123,43 +123,10 @@ set([hXLabel, hYLabel], 'FontSize', 14) set(hTitle, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title colorbar; -%% Heat map of angular spectral weight - -N_x = length(options.scan_groups); -N_y = length(results_all); - -BFields = [2.35, 2.15, 2.0, 1.85, 1.7, 1.55, 1.4, 1.35]; - -% Preallocate -angular_spectral_weight_matrix = zeros(N_y, N_x); - -for i = 1:N_y - for j = 1:N_x - angular_spectral_weight_matrix(i, j) = results_all(i).angular_spectral_weight(j); - end -end - -% Plot heatmap - -font = 'Bahnschrift'; - -figure(2) -clf -set(gcf,'Position',[50 50 950 750]) -imagesc(options.scan_groups, BFields, angular_spectral_weight_matrix); -colormap(sky); -set(gca, 'FontSize', 14, 'YDir', 'normal'); -hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); -hYLabel = ylabel('BField (G)', 'Interpreter', 'tex'); -hTitle = title('Angular Spectral Weight'); -set([hXLabel, hYLabel], 'FontSize', 14) -set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title -colorbar; - %% Heat map of radial spectral contrast -N_x = length(options.scan_groups); -N_y = length(results_all); +N_x = length(options.scan_groups); +N_y = length(results_all); BFields = [2.35, 2.15, 2.0, 1.85, 1.7, 1.55, 1.4, 1.35]; @@ -174,7 +141,7 @@ end % Plot heatmap -font = 'Bahnschrift'; +font = 'Bahnschrift'; figure(3) clf