From 96c0a50beba8d113a9bb1f75aa74e95f19321a96 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Fri, 11 Jul 2025 16:33:03 +0200 Subject: [PATCH] Changes to analyze datasets --- Data-Analyzer/analyzeFolder.m | 64 +++++++++++++------------ Data-Analyzer/conductSpectralAnalysis.m | 10 ++-- Data-Analyzer/extractQuantities.m | 14 +++--- 3 files changed, 46 insertions(+), 42 deletions(-) diff --git a/Data-Analyzer/analyzeFolder.m b/Data-Analyzer/analyzeFolder.m index cf4b1ca..a6a980e 100644 --- a/Data-Analyzer/analyzeFolder.m +++ b/Data-Analyzer/analyzeFolder.m @@ -180,7 +180,11 @@ function results = analyzeFolder(options) for k = 1:N_shots IMG = od_imgs{k}; - [IMGFFT, ~] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization); + if ~(max(IMG(:)) > 1) + IMGFFT = NaN(size(IMG)); + else + [IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization); + end % Size of original image (in pixels) [Ny, Nx] = size(IMG); @@ -214,7 +218,7 @@ function results = analyzeFolder(options) 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, []); + [theta_vals, S_theta] = computeAngularSpectralDistribution(fft_imgs{k}, kx, ky, k_min, k_max, N_angular_bins, Angular_Threshold, Angular_Sigma, []); [k_rho_vals, S_k] = computeRadialSpectralDistribution(fft_imgs{k}, kx, ky, theta_min, theta_max, N_radial_bins); S_k_smoothed = movmean(S_k, Radial_WindowSize); % Compute moving average (use convolution) or use conv for more control spectral_distribution{k} = S_theta; @@ -234,8 +238,8 @@ function results = analyzeFolder(options) % Loop through each unique theta and compute mean and standard error for i = 1:length(unique_scan_parameter_values) group_vals = radial_spectral_contrast(idx == i); - mean_rsc(i) = mean(group_vals); - stderr_rsc(i) = std(group_vals) / sqrt(length(group_vals)); % standard error = std / sqrt(N) + mean_rsc(i) = mean(group_vals, 'omitnan'); + stderr_rsc(i) = std(group_vals, 'omitnan') / sqrt(length(group_vals)); % standard error = std / sqrt(N) end % Preallocate arrays @@ -245,8 +249,8 @@ function results = analyzeFolder(options) % Loop through each unique theta and compute mean and standard error for i = 1:length(unique_scan_parameter_values) group_vals = angular_spectral_weight(idx == i); - mean_asw(i) = mean(group_vals); - stderr_asw(i) = std(group_vals) / sqrt(length(group_vals)); % standard error = std / sqrt(N) + mean_asw(i) = mean(group_vals, 'omitnan'); + stderr_asw(i) = std(group_vals, 'omitnan') / sqrt(length(group_vals)); % standard error = std / sqrt(N) end % Convert spectral distribution to matrix (N_shots x N_angular_bins) @@ -314,8 +318,8 @@ function results = analyzeFolder(options) 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); + num = mean(ref .* sec, 'omitnan'); + denom = mean(ref.^2, 'omitnan'); g2 = num / denom; correlations(k) = g2; @@ -456,52 +460,50 @@ function [k_rho_vals, S_radial] = computeRadialSpectralDistribution(IMGFFT, kx, end end -function [theta_vals, S_theta] = computeAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma, windowSize) +function [theta_vals, S_theta] = computeAngularSpectralDistribution(IMGFFT, kx, ky, k_min, k_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] + % Create wavenumber meshgrid + [KX, KY] = meshgrid(kx, ky); + Kmag = sqrt(KX.^2 + KY.^2); % radial wavenumber magnitude + Theta = atan2(KY, KX); % range [-pi, pi] - % Choose radial band - radial_mask = (R >= r_min) & (R <= r_max); + % Restrict to radial band in wavenumber space + radial_mask = (Kmag >= k_min) & (Kmag <= k_max); % Initialize angular structure factor - S_theta = zeros(1, num_bins); - theta_vals = linspace(0, pi, num_bins); + S_theta = zeros(1, num_bins); + theta_vals = linspace(0, pi, num_bins); % only 0 to pi due to symmetry - % Loop through angle bins + % Loop over angular bins for i = 1:num_bins - angle_start = (i-1) * pi / num_bins; + angle_start = (i - 1) * pi / num_bins; angle_end = i * pi / num_bins; - angle_mask = (Theta >= angle_start & Theta < angle_end); + 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 + % Optional smoothing if exist('sigma', 'var') && ~isempty(sigma) - % Gaussian convolution + % Gaussian smoothing 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)], ... + 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); + 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); + % Moving average smoothing + 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); + 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 diff --git a/Data-Analyzer/conductSpectralAnalysis.m b/Data-Analyzer/conductSpectralAnalysis.m index d31ff8c..839cf77 100644 --- a/Data-Analyzer/conductSpectralAnalysis.m +++ b/Data-Analyzer/conductSpectralAnalysis.m @@ -179,7 +179,7 @@ else end end -% ===== Get rotation angles ===== +%% ===== Get rotation angles ===== scan_parameter_values = zeros(1, length(files)); % Get information about the '/globals' group @@ -264,8 +264,12 @@ s_theta_list = cell(1, N_shots); % Angular spectrum for k = 1:N_shots IMG = od_imgs{k}; - [IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization); - + if ~(max(IMG(:)) > 1) + IMGFFT = NaN(size(IMG)); + else + [IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization); + end + % Size of original image (in pixels) [Ny, Nx] = size(IMG); diff --git a/Data-Analyzer/extractQuantities.m b/Data-Analyzer/extractQuantities.m index f75a9b7..dccb69f 100644 --- a/Data-Analyzer/extractQuantities.m +++ b/Data-Analyzer/extractQuantities.m @@ -2,8 +2,6 @@ % === Define folders and settings === -% === Define folders and settings === - baseFolder = '//DyLabNAS/Data/TwoDGas/2025/04/'; dates = ["01", "02"]; % Example: three folders @@ -39,8 +37,8 @@ options.Radial_WindowSize = 5; % Must be odd % Fourier analysis: Angular options.r_min = 10; options.r_max = 20; -options.k_min = 1.2; % in μm⁻¹ -options.k_max = 2.2; % in μm⁻¹ +options.k_min = 1.2771; % in μm⁻¹ +options.k_max = 2.5541; % in μm⁻¹ options.N_angular_bins = 180; options.Angular_Threshold = 75; options.Angular_Sigma = 2; @@ -125,7 +123,7 @@ set([hXLabel, hYLabel], 'FontSize', 14) set(hTitle, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title colorbar; -%% HEat map of angular spectral weight +%% Heat map of angular spectral weight N_x = length(options.scan_groups); N_y = length(results_all); @@ -158,7 +156,7 @@ 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 +%% Heat map of radial spectral contrast N_x = length(options.scan_groups); N_y = length(results_all); @@ -178,12 +176,12 @@ end font = 'Bahnschrift'; -figure(2) +figure(3) clf set(gcf,'Position',[50 50 950 750]) imagesc(options.scan_groups, BFields, radial_spectral_contrast_matrix); colormap(sky); -clim([0 0.005]) +clim([0 0.008]) set(gca, 'FontSize', 14, 'YDir', 'normal'); hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); hYLabel = ylabel('BField (G)', 'Interpreter', 'tex');