clear all close all %% 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 = "E:/Data - Experiment/2025/07/04/"; run = '0016'; folderPath = strcat(folderPath, run); cam = 5; angle = 0; center = [1430, 2040]; span = [200, 200]; fraction = [0.1, 0.1]; pixel_size = 5.86e-6; removeFringes = false; % 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 % scan_parameter = 'ps_rot_mag_fin_pol_angle'; scan_parameter = 'rot_mag_field'; % scan_parameter_text = 'Angle = '; scan_parameter_text = 'BField = '; savefolderPath = 'E:/Results - Experiment/B2.35G/'; savefileName = 'Droplets'; font = 'Bahnschrift'; skipUnshuffling = true; if strcmp(savefileName, 'DropletsToStripes') scan_groups = 0:5:45; elseif strcmp(savefileName, 'StripesToDroplets') scan_groups = 45:-5:0; end 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. 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; fullFileName = fullfile(files(k).folder, baseFileName); fprintf(1, 'Now reading %s\n', fullFileName); 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)); refimages(:,:,k) = subtractBackgroundOffset(cropODImage(bkg_img, center, span), fraction)'; absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img), center, span), fraction)'; end % 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 %% Get rotation angles scan_parameter_values = zeros(1, length(files)); % 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 %% Unshuffle if necessary to do so 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 %% Run Fourier analysis over images fft_imgs = cell(1, nimgs); spectral_contrast = zeros(1, nimgs); spectral_weight = zeros(1, nimgs); N_shots = length(od_imgs); % Create VideoWriter object for movie videoFile = VideoWriter([savefileName '.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:N_shots IMG = od_imgs{k}; [IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization); [rows, cols] = size(IMGFFT); mid_x = floor(cols/2); mid_y = floor(rows/2); fft_imgs{k} = IMGFFT(mid_y-zoom_size:mid_y+zoom_size, mid_x-zoom_size:mid_x+zoom_size); [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(fft_imgs{k}, r_min, r_max, N_angular_bins, Angular_Threshold, Angular_Sigma, []); [k_rho_vals, S_k] = computeRadialSpectralDistribution(fft_imgs{k}, 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_contrast(k) = computeSpectralContrast(fft_imgs{k}, r_min, r_max, Angular_Threshold); spectral_weight(k) = trapz(theta_vals, S_theta); figure(1); clf set(gcf,'Position',[500 100 1000 800]) t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid % Calculate the x and y limits for the cropped image y_min = center(1) - span(2) / 2; y_max = center(1) + span(2) / 2; x_min = center(2) - span(1) / 2; x_max = center(2) + span(1) / 2; % Generate x and y arrays representing the original coordinates for each pixel x_range = linspace(x_min, x_max, span(1)); y_range = linspace(y_min, y_max, span(2)); % 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, [scan_parameter_text, num2str(scan_parameter_values(k), '%.2f')], ... 'Color', 'white', 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontSize', 20, 'Units', 'normalized', 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top'); axis equal tight; hcb = colorbar; colormap(ax1, 'hot'); 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 (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) set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold'); % Set font and size for title % Plot the power spectrum ax2 = nexttile; imagesc(log(1 + abs(fft_imgs{k}).^2)); % Compute center of the FFT image [ny, nx] = size(fft_imgs{k}); cx = ceil(nx/2); cy = ceil(ny/2); % Define angles for the circle theta = linspace(0, 2*pi, 500); % Circle 1 at r_min x1 = cx + r_min * cos(theta); y1 = cy + r_min * sin(theta); % Circle 2 at r_max x2 = cx + r_max * cos(theta); y2 = cy + r_max * sin(theta); % Plot the circles hold on; plot(x1, y1, 'w--', 'LineWidth', 1.0); % Cyan for r_min plot(x2, y2, 'w--', 'LineWidth', 1.0); % Magenta for r_max plot([1, nx], [cy, cy], 'w--', 'LineWidth', 1.0); % white dashed horizontal line hold off; % 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 equal tight; hcb = colorbar; colormap(ax2, Colormaps.inferno()); set(gca, 'FontSize', 14); % For tick labels only 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)', '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 smoothed radial distribution nexttile plot(k_rho_vals, S_k_smoothed, 'LineWidth', 2); set(gca, 'FontSize', 14); % Tick labels set(gca, 'YScale', 'log'); % Logarithmic y-axis xlim([min(k_rho_vals), max(k_rho_vals)]); ylim([1, 1E8]) hXLabel = xlabel('k_\rho', 'Interpreter', 'tex'); hYLabel = ylabel('Magnitude (a.u.)', 'Interpreter', 'tex'); hTitle = title('Radial Spectral Distribution - S(k)', 'Interpreter', 'tex'); set([hXLabel, hYLabel], 'FontSize', 14, 'FontName', font); set(hTitle, 'FontSize', 16, 'FontWeight', 'bold', 'FontName', font); grid on; % Plot the angular distribution nexttile plot(theta_vals/pi, S_theta,'Linewidth',2); set(gca, 'FontSize', 14); % For tick labels only hXLabel = xlabel('\theta/\pi [rad]', 'Interpreter', 'tex'); hYLabel = ylabel('Normalized magnitude (a.u.)', '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 grid on 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); %% Helper Functions 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. % % Inputs: % I - Grayscale or RGB image matrix % % 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 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 [k_rho_vals, S_radial] = computeRadialSpectralDistribution(IMGFFT, thetamin, thetamax, num_bins) % IMGFFT : 2D FFT (should be fftshifted already) % thetamin : Minimum angle (in radians) % thetamax : Maximum angle (in radians) % num_radial_bins : Number of radial bins % sigma : Gaussian smoothing width (in bins) % Image size and center [ny, nx] = size(IMGFFT); [X, Y] = meshgrid(1:nx, 1:ny); cx = ceil(nx / 2); cy = ceil(ny / 2); dX = X - cx; dY = Y - cy; % Polar coordinates R = sqrt(dX.^2 + dY.^2); % radial coordinate Theta = atan2(dY, dX); % angle in radians [-pi, pi] % Angular mask (support wraparound from +pi to -pi) if thetamin < thetamax angle_mask = (Theta >= thetamin) & (Theta <= thetamax); else angle_mask = (Theta >= thetamin) | (Theta <= thetamax); end % Define full radial range: from center to farthest corner r_min = 0; r_max = sqrt((max([cx-1, nx-cx]))^2 + (max([cy-1, ny-cy]))^2); % Radial bins r_edges = linspace(r_min, r_max, num_bins + 1); k_rho_vals = 0.5 * (r_edges(1:end-1) + r_edges(2:end)); S_radial = zeros(1, num_bins); % Power spectrum power_spectrum = abs(IMGFFT).^2; % Radial integration over selected angles for i = 1:num_bins r_low = r_edges(i); r_high = r_edges(i + 1); radial_mask = (R >= r_low) & (R < r_high); full_mask = radial_mask & angle_mask; S_radial(i) = sum(power_spectrum(full_mask)); end end function [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, r_min, r_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] % Choose radial band radial_mask = (R >= r_min) & (R <= r_max); % Initialize angular structure factor S_theta = zeros(1, num_bins); theta_vals = linspace(0, pi, num_bins); % Loop through angle bins for i = 1:num_bins angle_start = (i-1) * pi / num_bins; angle_end = i * pi / num_bins; 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 if exist('sigma', 'var') && ~isempty(sigma) % Gaussian convolution 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)], ... gauss_kernel, 'same'); 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); 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); end % Normalize S_theta = S_theta / max(S_theta); end function contrast = computeSpectralContrast(IMGFFT, r_min, r_max, threshold) % 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); % Ring region (annulus) mask ring_mask = (R >= r_min) & (R <= r_max); % Squared magnitude in the ring ring_power = abs(IMGFFT).^2 .* ring_mask; % Maximum power in the ring ring_max = max(ring_power(:)); % Power at the DC component dc_power = abs(IMGFFT(cy, cx))^2; % Avoid division by zero if dc_power == 0 contrast = Inf; % or NaN or 0, depending on how you want to handle this else contrast = ring_max / dc_power; end end 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 ret = calculateODImage(imageAtom, imageBackground, imageDark) % Calculate the OD image for absorption imaging. % :param imageAtom: The image with atoms % :type imageAtom: numpy array % :param imageBackground: The image without atoms % :type imageBackground: numpy array % :param imageDark: The image without light % :type imageDark: numpy array % :return: The OD images % :rtype: numpy array numerator = imageBackground - imageDark; denominator = imageAtom - imageDark; numerator(numerator == 0) = 1; denominator(denominator == 0) = 1; ret = -log(double(abs(denominator ./ numerator))); if numel(ret) == 1 ret = ret(1); 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