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

Minor modifications to labels, changed computation of spectral weight.

Browse Source
This commit is contained in:
Karthik 2025-04-11 14:58:50 +02:00
parent d4c0dc3f07
commit 88febbd045
1 changed files with 44 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

89
Data-Analyzer/fourierAnalysis.m
View File

@ -1,14 +1,14 @@
%% 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/";
run = '0013';
folderPath = "C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/"; folderPath = strcat(folderPath, run);
run = '0013';
folderPath = strcat(folderPath, run);
cam = 5; cam = 5;
@ -23,10 +23,10 @@ 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.
filePattern = fullfile(folderPath, '*.h5'); filePattern = fullfile(folderPath, '*.h5');
files = dir(filePattern); files = dir(filePattern);
refimages = zeros(span(1) + 1, span(2) + 1, length(files)); refimages = zeros(span(1) + 1, span(2) + 1, length(files));
absimages = zeros(span(1) + 1, span(2) + 1, length(files)); absimages = zeros(span(1) + 1, span(2) + 1, length(files));
for k = 1 : length(files) for k = 1 : length(files)
baseFileName = files(k).name; baseFileName = files(k).name;
@ -80,7 +80,7 @@ end
%% Run Fourier analysis over images %% Run Fourier analysis over images
fft_imgs = cell(1, nimgs); fft_imgs = cell(1, nimgs);
integrated_sf = zeros(1, nimgs); spectral_weight = zeros(1, nimgs);
% Create VideoWriter object for movie % Create VideoWriter object for movie
videoFile = VideoWriter('Single_Shot_FFT.mp4', 'MPEG-4'); videoFile = VideoWriter('Single_Shot_FFT.mp4', 'MPEG-4');
@ -127,8 +127,8 @@ for k = 1 : length(od_imgs)
set(gca,'YDir','normal') set(gca,'YDir','normal')
set(gca, 'YTick', linspace(y_min, y_max, 5)); % Define y ticks 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 set(gca, 'YTickLabel', flip(linspace(y_min, y_max, 5))); % Flip only the labels
hXLabel = xlabel('X (pixels)', 'Interpreter', 'tex'); hXLabel = xlabel('x (pixels)', 'Interpreter', 'tex');
hYLabel = ylabel('Y (pixels)', 'Interpreter', 'tex'); hYLabel = ylabel('y (pixels)', 'Interpreter', 'tex');
hTitle = title('OD Image', 'Interpreter', 'tex'); hTitle = title('OD Image', 'Interpreter', 'tex');
set([hXLabel, hYLabel, hL, hText], 'FontName', font) set([hXLabel, hYLabel, hL, hText], 'FontName', font)
set([hXLabel, hYLabel, hL], 'FontSize', 14) set([hXLabel, hYLabel, hL], 'FontSize', 14)
@ -143,8 +143,8 @@ for k = 1 : length(od_imgs)
set(gca,'YDir','normal') set(gca,'YDir','normal')
set(gca, 'YTick', linspace(y_min, y_max, 5)); % Define y ticks 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 set(gca, 'YTickLabel', flip(linspace(y_min, y_max, 5))); % Flip only the labels
hXLabel = xlabel('X (pixels)', 'Interpreter', 'tex'); hXLabel = xlabel('x (pixels)', 'Interpreter', 'tex');
hYLabel = ylabel('Y (pixels)', 'Interpreter', 'tex'); hYLabel = ylabel('y (pixels)', 'Interpreter', 'tex');
hTitle = title('Denoised - Masked - Binarized', 'Interpreter', 'tex'); hTitle = title('Denoised - Masked - Binarized', 'Interpreter', 'tex');
set([hXLabel, hYLabel], 'FontName', font) set([hXLabel, hYLabel], 'FontName', font)
@ -170,9 +170,9 @@ for k = 1 : length(od_imgs)
colormap(ax3, 'jet'); colormap(ax3, 'jet');
set(gca, 'FontSize', 14); % For tick labels only set(gca, 'FontSize', 14); % For tick labels only
set(gca,'YDir','normal') set(gca,'YDir','normal')
hXLabel = xlabel('X (pixels)', 'Interpreter', 'tex'); hXLabel = xlabel('k_x', 'Interpreter', 'tex');
hYLabel = ylabel('Y (pixels)', 'Interpreter', 'tex'); hYLabel = ylabel('k_y', 'Interpreter', 'tex');
hTitle = title('Fourier Power Spectrum', 'Interpreter', 'tex'); hTitle = title('Power Spectrum - |S(k_x,k_y)|^2', 'Interpreter', 'tex');
set([hXLabel, hYLabel, hText], 'FontName', font) set([hXLabel, hYLabel, hText], 'FontName', font)
set([hXLabel, hYLabel], 'FontSize', 14) set([hXLabel, hYLabel], 'FontSize', 14)
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title
@ -190,13 +190,13 @@ for k = 1 : length(od_imgs)
%} %}
% Plot the angular structure factor % Plot the angular structure factor
nexttile nexttile
[theta_vals, S_theta] = computeAngularStructureFactor(zoomedIMGFFT, 10, 20, 180, 75, 2); [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(zoomedIMGFFT, 10, 20, 180, 75, 2);
integrated_sf(k) = trapz(theta_vals, S_theta); spectral_weight(k) = trapz(theta_vals, sqrt(S_theta));
plot(theta_vals/pi, S_theta,'Linewidth',2); plot(theta_vals/pi, S_theta,'Linewidth',2);
set(gca, 'FontSize', 14); % For tick labels only set(gca, 'FontSize', 14); % For tick labels only
hXLabel = xlabel('\theta (\pi)', 'Interpreter', 'tex'); hXLabel = xlabel('\theta (\pi)', 'Interpreter', 'tex');
hYLabel = ylabel('S(\theta)', 'Interpreter', 'tex'); hYLabel = ylabel('Normalized magnitude (a.u.)', 'Interpreter', 'tex');
hTitle = title('Angular Structure Factor', 'Interpreter', 'tex'); hTitle = title('Angular Spectral Distribution - |S(\theta)|^2', 'Interpreter', 'tex');
set([hXLabel, hYLabel, hText], 'FontName', font) set([hXLabel, hYLabel, hText], 'FontName', font)
set([hXLabel, hYLabel], 'FontSize', 14) set([hXLabel, hYLabel], 'FontSize', 14)
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title
@ -213,9 +213,9 @@ end
% Close the video file % Close the video file
close(videoFile); close(videoFile);
%% Track integrated structure factor across the transition %% Track spectral weight across the transition
% Assuming theta_values and integrated_sf are column vectors (or row vectors of same length) % Assuming theta_values and spectral_weight are column vectors (or row vectors of same length)
[unique_theta, ~, idx] = unique(theta_values); [unique_theta, ~, idx] = unique(theta_values);
% Preallocate arrays % Preallocate arrays
@ -224,7 +224,7 @@ stderr_sf = zeros(size(unique_theta));
% Loop through each unique theta and compute mean and standard error % Loop through each unique theta and compute mean and standard error
for i = 1:length(unique_theta) for i = 1:length(unique_theta)
group_vals = integrated_sf(idx == i); group_vals = spectral_weight(idx == i);
mean_sf(i) = mean(group_vals); mean_sf(i) = mean(group_vals);
stderr_sf(i) = std(group_vals) / sqrt(length(group_vals)); % standard error = std / sqrt(N) stderr_sf(i) = std(group_vals) / sqrt(length(group_vals)); % standard error = std / sqrt(N)
end end
@ -235,7 +235,7 @@ errorbar(unique_theta, mean_sf, stderr_sf, 'o--', ...
'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5); 'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5);
set(gca, 'FontSize', 14); % For tick labels only set(gca, 'FontSize', 14); % For tick labels only
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
hYLabel = ylabel('Integrated S(\theta)', 'Interpreter', 'tex'); hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex');
hTitle = title('Change during rotation', 'Interpreter', 'tex'); hTitle = title('Change during rotation', 'Interpreter', 'tex');
set([hXLabel, hYLabel], 'FontName', font) set([hXLabel, hYLabel], 'FontName', font)
set([hXLabel, hYLabel], 'FontSize', 14) set([hXLabel, hYLabel], 'FontSize', 14)
@ -245,7 +245,7 @@ grid on
%% k-means Clustering %% k-means Clustering
% Reshape to column vector % Reshape to column vector
X = mean_sf(:); X = mean_sf(:);
% Determine the number of clusters to try (you can experiment with different values here) % Determine the number of clusters to try (you can experiment with different values here)
optimalClusters = 3; % Based on prior knowledge or experimentation optimalClusters = 3; % Based on prior knowledge or experimentation
@ -254,10 +254,10 @@ optimalClusters = 3; % Based on prior knowledge or experimentation
rng(42); rng(42);
% Specify initialization method ('plus' is the default) % Specify initialization method ('plus' is the default)
startMethod = 'plus'; % Options: 'uniform', 'plus', 'sample' startMethod = 'plus'; % Options: 'uniform', 'plus', 'sample'
% Apply K-means clustering with controlled initialization % Apply K-means clustering with controlled initialization
[idx, C] = kmeans(X, optimalClusters, 'Start', startMethod); [idx, C] = kmeans(X, optimalClusters, 'Start', startMethod);
% Plot the results % Plot the results
figure(3); figure(3);
@ -272,10 +272,10 @@ errorbar(unique_theta, mean_sf, stderr_sf, 'o--', ...
scatter(unique_theta, X, 50, idx, 'filled'); scatter(unique_theta, X, 50, idx, 'filled');
% Get the current y-axis limits % Get the current y-axis limits
current_ylim = ylim; current_ylim = ylim;
% Generate colors for each cluster % Generate colors for each cluster
colors = lines(optimalClusters); % Create distinct colors for each cluster colors = lines(optimalClusters); % Create distinct colors for each cluster
% Loop through each cluster and fill the regions % Loop through each cluster and fill the regions
for i = 1:optimalClusters for i = 1:optimalClusters
@ -293,7 +293,7 @@ for i = 1:optimalClusters
end end
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
hYLabel = ylabel('Integrated S(\theta)', 'Interpreter', 'tex'); hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex');
hTitle = title('Regimes identified by k-means clustering', 'Interpreter', 'tex'); hTitle = title('Regimes identified by k-means clustering', 'Interpreter', 'tex');
set([hXLabel, hYLabel], 'FontName', font) set([hXLabel, hYLabel], 'FontName', font)
set([hXLabel, hYLabel], 'FontSize', 14) set([hXLabel, hYLabel], 'FontSize', 14)
@ -379,7 +379,7 @@ function [IMGFFT, IMGBIN] = computeFourierTransform(I)
end end
function [theta_vals, S_theta] = computeAngularStructureFactor(IMGFFT, r_min, r_max, num_bins, threshold, sigma) function [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma)
% Apply threshold to isolate strong peaks % Apply threshold to isolate strong peaks
IMGFFT(IMGFFT < threshold) = 0; IMGFFT(IMGFFT < threshold) = 0;
@ -402,33 +402,32 @@ function [theta_vals, S_theta] = computeAngularStructureFactor(IMGFFT, r_min, r_
% Loop through each angle bin % Loop through each angle bin
for i = 1:180 for i = 1:180
angle_start = (i-1) * pi / num_bins; angle_start = (i-1) * pi / num_bins;
angle_end = i * pi / num_bins; angle_end = i * pi / num_bins;
% Define a mask for the given angle range % Define a mask for the given angle range
angle_mask = (Theta >= angle_start & Theta < angle_end); angle_mask = (Theta >= angle_start & Theta < angle_end);
bin_mask = radial_mask & angle_mask; bin_mask = radial_mask & angle_mask;
% Extract the Fourier components for the given angle % Extract the Fourier components for the given angle
fft_angle = IMGFFT .* bin_mask; fft_angle = IMGFFT .* bin_mask;
% Integrate the Fourier components over the radius at the angle % Integrate the Fourier components over the radius at the angle
S_theta(i) = sum(sum(abs(fft_angle).^2)); % Compute structure factor (sum of squared magnitudes) S_theta(i) = sum(sum(abs(fft_angle).^2)); % sum of squared magnitudes
end end
% Create a 1D Gaussian kernel % Create a 1D Gaussian kernel
half_width = ceil(3 * sigma); half_width = ceil(3 * sigma);
x = -half_width:half_width; x = -half_width:half_width;
gauss_kernel = exp(-x.^2 / (2 * sigma^2)); gauss_kernel = exp(-x.^2 / (2 * sigma^2));
gauss_kernel = gauss_kernel / sum(gauss_kernel); % normalize gauss_kernel = gauss_kernel / sum(gauss_kernel); % normalize
% Apply convolution (circular padding to preserve periodicity) % Apply convolution (circular padding to preserve periodicity)
S_theta = conv([S_theta(end-half_width+1:end), S_theta, S_theta(1:half_width)], gauss_kernel, 'same'); 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); % crop back to original size S_theta = S_theta(half_width+1:end-half_width); % crop back to original size
% Normalize to maximum value of 1 % Normalize to 1
S_theta = S_theta / max(S_theta); S_theta = S_theta / max(S_theta);
end end
function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction) function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction)