diff --git a/Data-Analyzer/fourierAnalysis.m b/Data-Analyzer/fourierAnalysis.m index a29b3d3..fc5aefc 100644 --- a/Data-Analyzer/fourierAnalysis.m +++ b/Data-Analyzer/fourierAnalysis.m @@ -1,14 +1,14 @@ %% 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/"; - -run = '0013'; - -folderPath = strcat(folderPath, run); +folderPath = strcat(folderPath, run); cam = 5; @@ -23,10 +23,10 @@ removeFringes = false; %% 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)); +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; @@ -80,7 +80,7 @@ end %% Run Fourier analysis over images fft_imgs = cell(1, nimgs); -integrated_sf = zeros(1, nimgs); +spectral_weight = zeros(1, nimgs); % Create VideoWriter object for movie videoFile = VideoWriter('Single_Shot_FFT.mp4', 'MPEG-4'); @@ -127,8 +127,8 @@ for k = 1 : length(od_imgs) 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 (pixels)', 'Interpreter', 'tex'); - hYLabel = ylabel('Y (pixels)', 'Interpreter', 'tex'); + hXLabel = xlabel('x (pixels)', 'Interpreter', 'tex'); + hYLabel = ylabel('y (pixels)', 'Interpreter', 'tex'); hTitle = title('OD Image', 'Interpreter', 'tex'); set([hXLabel, hYLabel, hL, hText], 'FontName', font) set([hXLabel, hYLabel, hL], 'FontSize', 14) @@ -143,8 +143,8 @@ for k = 1 : length(od_imgs) 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 (pixels)', 'Interpreter', 'tex'); - hYLabel = ylabel('Y (pixels)', 'Interpreter', 'tex'); + hXLabel = xlabel('x (pixels)', 'Interpreter', 'tex'); + hYLabel = ylabel('y (pixels)', 'Interpreter', 'tex'); hTitle = title('Denoised - Masked - Binarized', 'Interpreter', 'tex'); set([hXLabel, hYLabel], 'FontName', font) @@ -170,9 +170,9 @@ for k = 1 : length(od_imgs) colormap(ax3, 'jet'); set(gca, 'FontSize', 14); % For tick labels only set(gca,'YDir','normal') - hXLabel = xlabel('X (pixels)', 'Interpreter', 'tex'); - hYLabel = ylabel('Y (pixels)', 'Interpreter', 'tex'); - hTitle = title('Fourier Power Spectrum', 'Interpreter', 'tex'); + hXLabel = xlabel('k_x', 'Interpreter', 'tex'); + hYLabel = ylabel('k_y', 'Interpreter', 'tex'); + hTitle = title('Power Spectrum - |S(k_x,k_y)|^2', 'Interpreter', 'tex'); set([hXLabel, hYLabel, hText], 'FontName', font) set([hXLabel, hYLabel], 'FontSize', 14) 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 nexttile - [theta_vals, S_theta] = computeAngularStructureFactor(zoomedIMGFFT, 10, 20, 180, 75, 2); - integrated_sf(k) = trapz(theta_vals, S_theta); + [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(zoomedIMGFFT, 10, 20, 180, 75, 2); + spectral_weight(k) = trapz(theta_vals, sqrt(S_theta)); plot(theta_vals/pi, S_theta,'Linewidth',2); set(gca, 'FontSize', 14); % For tick labels only hXLabel = xlabel('\theta (\pi)', 'Interpreter', 'tex'); - hYLabel = ylabel('S(\theta)', 'Interpreter', 'tex'); - hTitle = title('Angular Structure Factor', 'Interpreter', 'tex'); + hYLabel = ylabel('Normalized magnitude (a.u.)', 'Interpreter', 'tex'); + hTitle = title('Angular Spectral Distribution - |S(\theta)|^2', 'Interpreter', 'tex'); set([hXLabel, hYLabel, hText], 'FontName', font) set([hXLabel, hYLabel], 'FontSize', 14) set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title @@ -213,9 +213,9 @@ end % Close the video file 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); % Preallocate arrays @@ -224,7 +224,7 @@ stderr_sf = zeros(size(unique_theta)); % Loop through each unique theta and compute mean and standard error for i = 1:length(unique_theta) - group_vals = integrated_sf(idx == i); + group_vals = spectral_weight(idx == i); mean_sf(i) = mean(group_vals); stderr_sf(i) = std(group_vals) / sqrt(length(group_vals)); % standard error = std / sqrt(N) end @@ -235,7 +235,7 @@ errorbar(unique_theta, mean_sf, stderr_sf, 'o--', ... 'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5); set(gca, 'FontSize', 14); % For tick labels only 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'); set([hXLabel, hYLabel], 'FontName', font) set([hXLabel, hYLabel], 'FontSize', 14) @@ -245,7 +245,7 @@ grid on %% k-means Clustering % Reshape to column vector -X = mean_sf(:); +X = mean_sf(:); % Determine the number of clusters to try (you can experiment with different values here) optimalClusters = 3; % Based on prior knowledge or experimentation @@ -254,10 +254,10 @@ optimalClusters = 3; % Based on prior knowledge or experimentation rng(42); % 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 -[idx, C] = kmeans(X, optimalClusters, 'Start', startMethod); +[idx, C] = kmeans(X, optimalClusters, 'Start', startMethod); % Plot the results figure(3); @@ -272,10 +272,10 @@ errorbar(unique_theta, mean_sf, stderr_sf, 'o--', ... scatter(unique_theta, X, 50, idx, 'filled'); % Get the current y-axis limits -current_ylim = ylim; +current_ylim = ylim; % 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 for i = 1:optimalClusters @@ -293,7 +293,7 @@ for i = 1:optimalClusters end 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'); set([hXLabel, hYLabel], 'FontName', font) set([hXLabel, hYLabel], 'FontSize', 14) @@ -379,7 +379,7 @@ function [IMGFFT, IMGBIN] = computeFourierTransform(I) 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 IMGFFT(IMGFFT < threshold) = 0; @@ -402,33 +402,32 @@ function [theta_vals, S_theta] = computeAngularStructureFactor(IMGFFT, r_min, r_ % Loop through each angle bin for i = 1:180 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 - 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 - fft_angle = IMGFFT .* bin_mask; + fft_angle = IMGFFT .* bin_mask; % 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 % Create a 1D Gaussian kernel - half_width = ceil(3 * sigma); - x = -half_width:half_width; + 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); % normalize % 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 = S_theta(half_width+1:end-half_width); % crop back to original size - - % Normalize to maximum value of 1 + 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 + + % Normalize to 1 S_theta = S_theta / max(S_theta); - end function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction)