From 4918ee1ec079b961f1f30a7e074b1dad75cd92f4 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Wed, 23 Apr 2025 15:33:57 +0200 Subject: [PATCH] Updated Fourier analysis script, image plotting script. --- Data-Analyzer/fourierAnalysis.m | 291 ++++++++++++-------------------- Data-Analyzer/plotImages.m | 8 +- 2 files changed, 114 insertions(+), 185 deletions(-) diff --git a/Data-Analyzer/fourierAnalysis.m b/Data-Analyzer/fourierAnalysis.m index fc5aefc..ab26dd7 100644 --- a/Data-Analyzer/fourierAnalysis.m +++ b/Data-Analyzer/fourierAnalysis.m @@ -4,22 +4,35 @@ groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axi "/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/"; +folderPath = "C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/15042025/"; -run = '0013'; +run = '0035'; folderPath = strcat(folderPath, run); cam = 5; angle = 0; -center = [1285, 2105]; +center = [1300, 2108]; span = [200, 200]; fraction = [0.1, 0.1]; pixel_size = 5.86e-6; removeFringes = false; +% scan_parameter = 'rot_mag_fin_pol_angle'; +% scan_parameter = 'rot_mag_field'; +scan_parameter = 'rot_mag_field_up'; +% scan_parameter_text = 'Angle = '; +scan_parameter_text = 'BField = '; + +font = 'Bahnschrift'; + +skipPreprocessing = true; +skipMasking = true; +skipIntensityThresholding = true; +skipBinarization = true; + %% Compute OD image, rotate and extract ROI for analysis % Get a list of all files in the folder with the desired file name pattern. @@ -63,7 +76,7 @@ else end %% Get rotation angles -theta_values = zeros(1, length(files)); +scan_parameter_values = zeros(1, length(files)); % Get information about the '/globals' group for k = 1 : length(files) @@ -71,8 +84,12 @@ for k = 1 : length(files) fullFileName = fullfile(files(k).folder, baseFileName); info = h5info(fullFileName, '/globals'); for i = 1:length(info.Attributes) - if strcmp(info.Attributes(i).Name, 'rot_mag_fin_pol_angle') - theta_values(k) = 180 - info.Attributes(i).Value; + if strcmp(info.Attributes(i).Name, scan_parameter) + if strcmp(scan_parameter, '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 @@ -91,14 +108,13 @@ open(videoFile); % Open the video file to write % Display the cropped image for k = 1 : length(od_imgs) IMG = od_imgs{k}; - [IMGFFT, IMGBIN] = computeFourierTransform(IMG); + [IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization); figure(1); clf set(gcf,'Position',[500 100 1000 800]) t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid - font = 'Bahnschrift'; - + % Calculate the x and y limits for the cropped image y_min = center(1) - span(2) / 2; y_max = center(1) + span(2) / 2; @@ -109,14 +125,14 @@ for k = 1 : length(od_imgs) x_range = linspace(x_min, x_max, span(1)); y_range = linspace(y_min, y_max, span(2)); - % Display the cropped image + % Display the cropped OD image ax1 = nexttile; imagesc(x_range, y_range, IMG) % 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 % Top-right corner (normalized axis coordinates) - hText = text(1 - x_offset, 1 - y_offset, ['Angle: ', num2str(theta_values(k), '%.1f')], ... + hText = text(1 - x_offset, 1 - y_offset, [scan_parameter_text, num2str(scan_parameter_values(k), '%.1f')], ... 'Color', 'white', 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontSize', 20, 'Units', 'normalized', 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top'); axis equal tight; hcb = colorbar; @@ -134,8 +150,9 @@ for k = 1 : length(od_imgs) set([hXLabel, hYLabel, hL], 'FontSize', 14) set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title + % Plot the processed image ax2 = nexttile; - imagesc(x_range, y_range, IMGBIN) + imagesc(x_range, y_range, IMGPR) axis equal tight; hcb = colorbar; colormap(ax2, 'parula'); @@ -145,26 +162,23 @@ for k = 1 : length(od_imgs) 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'); - hTitle = title('Denoised - Masked - Binarized', 'Interpreter', 'tex'); + hTitle = title('Processed Image', 'Interpreter', 'tex'); set([hXLabel, hYLabel], 'FontName', font) set([hXLabel, hYLabel], 'FontSize', 14) set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title + % Plot the power spectrum ax3 = nexttile; [rows, cols] = size(IMGFFT); zoom_size = 50; % Zoomed-in region around center mid_x = floor(cols/2); mid_y = floor(rows/2); - zoomedIMGFFT = IMGFFT(mid_y-zoom_size:mid_y+zoom_size, mid_x-zoom_size:mid_x+zoom_size); - fft_imgs{k} = zoomedIMGFFT; - imagesc(log(1 + zoomedIMGFFT)); + fft_imgs{k} = IMGFFT(mid_y-zoom_size:mid_y+zoom_size, mid_x-zoom_size:mid_x+zoom_size); + imagesc(log(1 + abs(fft_imgs{k}).^2)); % 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 - % Top-right corner (normalized axis coordinates) - % hText = text(1 - x_offset, 1 - y_offset, ['Angle: ', num2str(theta_values(k), '%.1f')], ... - % 'Color', 'white', 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontSize', 20, 'Units', 'normalized', 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top'); axis equal tight; hcb = colorbar; colormap(ax3, 'jet'); @@ -172,31 +186,20 @@ for k = 1 : length(od_imgs) set(gca,'YDir','normal') hXLabel = xlabel('k_x', 'Interpreter', 'tex'); hYLabel = ylabel('k_y', 'Interpreter', 'tex'); - hTitle = title('Power Spectrum - |S(k_x,k_y)|^2', 'Interpreter', 'tex'); + hTitle = title('Power Spectrum - S(k_x,k_y)', '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 % Plot the angular distribution - %{ nexttile - [theta_vals, angular_intensity] = computeAngularDistribution(zoomedIMGFFT, 10, 20, 100, 75); - polarhistogram('BinEdges', theta_vals, 'BinCounts', angular_intensity, ... - 'FaceColor', [0.2 0.6 0.9], 'EdgeColor', 'k'); - set(gca, 'FontSize', 14); % For tick labels only - hTitle = title('Angular Distribution', 'Interpreter', 'tex'); - set(hTitle, 'FontName', font) - set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title - %} - % Plot the angular structure factor - nexttile - [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(zoomedIMGFFT, 10, 20, 180, 75, 2); - spectral_weight(k) = trapz(theta_vals, sqrt(S_theta)); + [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(fft_imgs{k}, 10, 20, 180, 75, 2); + spectral_weight(k) = trapz(theta_vals, S_theta); plot(theta_vals/pi, S_theta,'Linewidth',2); set(gca, 'FontSize', 14); % For tick labels only - hXLabel = xlabel('\theta (\pi)', 'Interpreter', 'tex'); + hXLabel = xlabel('\theta/\pi [rad]', 'Interpreter', 'tex'); hYLabel = ylabel('Normalized magnitude (a.u.)', 'Interpreter', 'tex'); - hTitle = title('Angular Spectral Distribution - |S(\theta)|^2', 'Interpreter', 'tex'); + hTitle = title('Angular Spectral Distribution - S(\theta)', '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 @@ -215,15 +218,15 @@ close(videoFile); %% Track spectral weight across the transition -% Assuming theta_values and spectral_weight are column vectors (or row vectors of same length) -[unique_theta, ~, idx] = unique(theta_values); +% Assuming scan_parameter_values and spectral_weight are column vectors (or row vectors of same length) +[unique_scan_parameter_values, ~, idx] = unique(scan_parameter_values); % Preallocate arrays -mean_sf = zeros(size(unique_theta)); -stderr_sf = zeros(size(unique_theta)); +mean_sf = zeros(size(unique_scan_parameter_values)); +stderr_sf = zeros(size(unique_scan_parameter_values)); % Loop through each unique theta and compute mean and standard error -for i = 1:length(unique_theta) +for i = 1:length(unique_scan_parameter_values) 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) @@ -231,12 +234,13 @@ end figure(2); set(gcf,'Position',[100 100 950 750]) -errorbar(unique_theta, mean_sf, stderr_sf, 'o--', ... +errorbar(unique_scan_parameter_values, 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'); +% hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex'); +hXLabel = xlabel('B_z (G)', 'Interpreter', 'tex'); hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex'); -hTitle = title('Change during rotation', 'Interpreter', 'tex'); +hTitle = title('Change across transition', 'Interpreter', 'tex'); set([hXLabel, hYLabel], 'FontName', font) set([hXLabel, hYLabel], 'FontSize', 14) set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title @@ -265,11 +269,11 @@ set(gcf,'Position',[100 100 950 750]) hold on; % Plot error bars with mean_sf and stderr_sf -errorbar(unique_theta, mean_sf, stderr_sf, 'o--', ... +errorbar(unique_scan_parameter_values, mean_sf, stderr_sf, 'o--', ... 'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5); % Scatter plot for data points (showing clusters) -scatter(unique_theta, X, 50, idx, 'filled'); +scatter(unique_scan_parameter_values, X, 50, idx, 'filled'); % Get the current y-axis limits current_ylim = ylim; @@ -283,8 +287,8 @@ for i = 1:optimalClusters clusterIdx = find(idx == i); % Find the range of x-values for this cluster - x_min = unique_theta(clusterIdx(1)); % Starting x-value for the cluster - x_max = unique_theta(clusterIdx(end)); % Ending x-value for the cluster + x_min = unique_scan_parameter_values(clusterIdx(1)); % Starting x-value for the cluster + x_max = unique_scan_parameter_values(clusterIdx(end)); % Ending x-value for the cluster % Fill the region corresponding to the cluster fill([x_min, x_max, x_max, x_min], ... @@ -302,7 +306,7 @@ grid on; hold off; %% Helper Functions -function [IMGFFT, IMGBIN] = computeFourierTransform(I) +function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization) % computeFourierSpectrum - Computes the 2D Fourier power spectrum % of binarized and enhanced lattice image features, with optional central mask. % @@ -311,72 +315,66 @@ function [IMGFFT, IMGBIN] = computeFourierTransform(I) % % Output: % F_mag - 2D Fourier power spectrum (shifted) + + if ~skipPreprocessing + % Preprocessing: Denoise + filtered = imgaussfilt(I, 10); + IMGPR = I - filtered; % adjust sigma as needed + else + IMGPR = I; + end + + if ~skipMasking + [rows, cols] = size(IMGPR); + [X, Y] = meshgrid(1:cols, 1:rows); + % Elliptical mask parameters + cx = cols / 2; + cy = rows / 2; + + % Shifted coordinates + x = X - cx; + y = Y - cy; + + % Ellipse semi-axes + rx = 0.4 * cols; + ry = 0.2 * rows; + + % Rotation angle in degrees -> radians + theta_deg = 30; % Adjust as needed + theta = deg2rad(theta_deg); + + % Rotated ellipse equation + cos_t = cos(theta); + sin_t = sin(theta); + + x_rot = (x * cos_t + y * sin_t); + y_rot = (-x * sin_t + y * cos_t); + + ellipseMask = (x_rot.^2) / rx^2 + (y_rot.^2) / ry^2 <= 1; + + % Apply cutout mask + IMGPR = IMGPR .* ellipseMask; + end - % Preprocessing: Denoise - filtered = imgaussfilt(I, 10); - I_filt = I - filtered; % adjust sigma as needed - - % Elliptical mask parameters - [rows, cols] = size(I_filt); - [X, Y] = meshgrid(1:cols, 1:rows); - cx = cols / 2; - cy = rows / 2; - - % Shifted coordinates - x = X - cx; - y = Y - cy; - - % Ellipse semi-axes - rx = 0.4 * cols; - ry = 0.2 * rows; - - % Rotation angle in degrees -> radians - theta_deg = 30; % Adjust as needed - theta = deg2rad(theta_deg); - - % Rotated ellipse equation - cos_t = cos(theta); - sin_t = sin(theta); - - x_rot = (x * cos_t + y * sin_t); - y_rot = (-x * sin_t + y * cos_t); - - ellipseMask = (x_rot.^2) / rx^2 + (y_rot.^2) / ry^2 <= 1; - - % Apply cutout mask - I_masked = I_filt .* ellipseMask; - - % Apply global intensity threshold mask - intensity_thresh = 0.20; - intensity_mask = I_masked > intensity_thresh; - I_masked = I_masked .* intensity_mask; - - % Adaptive binarization and cleanup - IMGBIN = imbinarize(I_masked, 'adaptive', 'Sensitivity', 0.0); - IMGBIN = imdilate(IMGBIN, strel('disk', 2)); - IMGBIN = imerode(IMGBIN, strel('disk', 1)); - IMGBIN = imfill(IMGBIN, 'holes'); - - % Compute 2D Fourier Transform - F = fft2(double(I)); - IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center - - % Define the radius for the circular region to exclude - region_radius = 4; % Adjust the radius as needed - - % Create a circular mask - [~, center_idx] = max(IMGFFT(:)); - [cx, cy] = ind2sub(size(IMGFFT), center_idx); - - % Equation for a circle (centered at cx, cy) - center_region = (X - cx).^2 + (Y - cy).^2 <= region_radius^2; - - % Define a scaling factor for the central region (e.g., reduce amplitude by 90%) - scaling_factor = 0.1; % Scale center region by 10% - - % Apply the scaling factor to the center region - IMGFFT(center_region) = IMGFFT(center_region) * scaling_factor; + if ~skipIntensityThresholding + % Apply global intensity threshold mask + intensity_thresh = 0.20; + intensity_mask = IMGPR > intensity_thresh; + IMGPR = IMGPR .* intensity_mask; + end + if ~skipBinarization + % Adaptive binarization and cleanup + IMGPR = imbinarize(IMGPR, 'adaptive', 'Sensitivity', 0.0); + IMGPR = imdilate(IMGPR, strel('disk', 2)); + IMGPR = imerode(IMGPR, strel('disk', 1)); + IMGPR = imfill(IMGPR, 'holes'); + F = fft2(double(IMGPR)); % Compute 2D Fourier Transform + IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center + else + F = fft2(double(IMGPR)); % Compute 2D Fourier Transform + IMGFFT = abs(fftshift(F))'; % Shift zero frequency to center + end end function [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma) @@ -610,76 +608,7 @@ for k = 1 : length(od_imgs) end %} -%% Averaged FFT -%{ -% Assuming od_imgs is a cell array of size 4*n -n = length(fft_imgs) / 4; % Calculate n -fft_imgs_avg = cell(1, n); % Initialize the new cell array to hold the averaged images - -for i = 1:n - % Take the 4 corresponding images from od_imgs - img1 = fft_imgs{4*i-3}; % 1st image in the group - img2 = fft_imgs{4*i-2}; % 2nd image in the group - img3 = fft_imgs{4*i-1}; % 3rd image in the group - img4 = fft_imgs{4*i}; % 4th image in the group - - % Compute the average of these 4 images - avg_img = (img1 + img2 + img3 + img4) / 4; - - % Store the averaged image in the new cell array - fft_imgs_avg{i} = avg_img; -end - -% Create VideoWriter object for movie -videoFile = VideoWriter('Averaged_FFT.mp4', 'MPEG-4'); -videoFile.Quality = 100; % Set quality to maximum (0–100) -videoFile.FrameRate = 2; % Set the frame rate (frames per second) -open(videoFile); % Open the video file to write - -% Display the cropped image -for k = 1 : length(fft_imgs_avg) - figure(3) - clf - set(gcf,'Position',[50 50 1500 550]) - set(gca,'FontSize',16,'Box','On','Linewidth',2); - t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x2 grid - - nexttile - imagesc(log(1 + fft_imgs_avg{k})); - % 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 - % Top-right corner (normalized axis coordinates) - text(1 - x_offset, 1 - y_offset, ['Angle: ', num2str(theta_values(k), '%.1f')], ... - 'Color', 'white', 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontSize', 20, 'Units', 'normalized', 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top'); - axis equal tight; - hcb = colorbar; - set(gca,'YDir','normal') - xlabel('X', 'Interpreter', 'tex','FontSize',16); - ylabel('Y', 'Interpreter', 'tex','FontSize',16); - title('Averaged Fourier Power Spectrum','FontSize',16); - - - % Plot the angular structure factor - nexttile - [theta_vals, angular_intensity] = computeAngularDistribution(fft_imgs_avg{k}, 10, 20, 100, 75); - polarhistogram('BinEdges', theta_vals, 'BinCounts', angular_intensity, ... - 'FaceColor', [0.2 0.6 0.9], 'EdgeColor', 'k'); - title('Angular Distribution'); - - drawnow - pause(0.5) - - % Capture the current frame and write it to the video - frame = getframe(gcf); % Capture the current figure as a frame - writeVideo(videoFile, frame); % Write the frame to the video -end - -% Close the video file -close(videoFile); - -%} %% Angular Distribution %{ diff --git a/Data-Analyzer/plotImages.m b/Data-Analyzer/plotImages.m index 92a9bc1..9e97009 100644 --- a/Data-Analyzer/plotImages.m +++ b/Data-Analyzer/plotImages.m @@ -6,15 +6,15 @@ groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis_Ca folderPath = "C:/Users/Karthik/Documents/GitRepositories/Calculations/Data-Analyzer/"; -run = '0013'; +run = '0060'; folderPath = strcat(folderPath, run); cam = 5; angle = 0; -center = [1285, 2105]; -span = [200, 200]; +center = [1630, 1700]; +span = [500, 500]; fraction = [0.1, 0.1]; pixel_size = 5.86e-6; @@ -104,7 +104,7 @@ for k = 1 : length(od_imgs) hL = ylabel(hcb, 'Optical Density', 'FontSize', 16); set(hL,'Rotation',-90); colormap jet; - set(gca,'CLim',[0 3.0]); + set(gca,'CLim',[0 0.4]); 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