diff --git a/Data-Analyzer/conductSpectralAnalysis.m b/Data-Analyzer/conductSpectralAnalysis.m index 3d44afb..cbd173a 100644 --- a/Data-Analyzer/conductSpectralAnalysis.m +++ b/Data-Analyzer/conductSpectralAnalysis.m @@ -6,7 +6,7 @@ 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 = "E:/Data - Experiment/2025/07/04/"; +folderPath = "D:/Data - Experiment/2025/07/04/"; run = '0016'; @@ -20,8 +20,12 @@ span = [200, 200]; fraction = [0.1, 0.1]; pixel_size = 5.86e-6; +magnification = 23.94; removeFringes = false; +ImagingMode = 'HighIntensity'; +PulseDuration = 5e-6; % in s + % Fourier analysis settings % Radial Spectral Distribution @@ -61,6 +65,8 @@ skipPreprocessing = true; skipMasking = true; skipIntensityThresholding = true; skipBinarization = true; +skipMovieRender = true; +skipSaveFigures = 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. @@ -81,8 +87,7 @@ for k = 1 : length(files) 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)'; - + absimages(:,:,k) = subtractBackgroundOffset(cropODImage(calculateODImage(atm_img, bkg_img, dark_img, ImagingMode, PulseDuration), center, span), fraction)'; end % Fringe removal @@ -167,25 +172,68 @@ 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 +avg_ps_accum = 0; +avg_S_k_accum = 0; +avg_S_theta_accum = 0; + +% Pre-allocate once sizes are known (after first run) +fft_size_known = false; + +if ~skipMovieRender + % 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 +end + +if ~skipSaveFigures + % Define folder for saving images + saveFolder = [savefileName '_SavedFigures']; + if ~exist(saveFolder, 'dir') + mkdir(saveFolder); + end +end % 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); + % Size of original image (in pixels) + [Ny, Nx] = size(IMG); - mid_x = floor(cols/2); - mid_y = floor(rows/2); + % Real-space pixel size in micrometers after magnification + dx = pixel_size / magnification; + dy = dx; % assuming square pixels + + % Real-space axes + x = ((1:Nx) - ceil(Nx/2)) * dx * 1E6; + y = ((1:Ny) - ceil(Ny/2)) * dy * 1E6; + + % Reciprocal space increments (frequency domain, μm⁻¹) + dvx = 1 / (Nx * dx); + dvy = 1 / (Ny * dy); + + % Frequency axes + vx = (-floor(Nx/2):ceil(Nx/2)-1) * dvx; + vy = (-floor(Ny/2):ceil(Ny/2)-1) * dvy; + + % Wavenumber axes + kx_full = 2 * pi * vx * 1E-6; % μm⁻¹ + ky_full = 2 * pi * vy * 1E-6; + + % Crop FFT image around center + mid_x = floor(Nx/2); + mid_y = floor(Ny/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); + + % Crop wavenumber axes to match fft_imgs{k} + 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, []); + [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_contrast(k) = computeSpectralContrast(fft_imgs{k}, r_min, r_max, Angular_Threshold); @@ -194,121 +242,157 @@ for k = 1:N_shots 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 + t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); + + % ======= OD IMAGE (real space) ======= 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 + imagesc(x, y, IMG) + 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 + + % Convert pixel grid to µm (already done: x and y axes) + % Draw ↘ diagonal (top-left to bottom-right) + drawODOverlays(x(1), y(1), x(end), y(end)); + + % Draw ↙ diagonal (top-right to bottom-left) + drawODOverlays(x(end), y(1), x(1), y(end)); + 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; + set(gca, 'FontSize', 14, 'YDir', 'normal') + colormap(ax1, 'sky'); hcb = colorbar; + ylabel(hcb, 'Optical Density', 'Rotation', -90, 'FontSize', 14, 'FontName', font); + + xlabel('x (\mum)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font); + ylabel('y (\mum)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font); + title('OD Image', 'FontSize', 16, 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontName', font); + + text(0.975, 0.975, [scan_parameter_text, num2str(scan_parameter_values(k), '%.2f'), ' G'], ... + 'Color', 'black', 'FontWeight', 'bold', 'FontSize', 14, ... + 'Interpreter', 'tex', 'Units', 'normalized', ... + 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top'); + + % ======= FFT POWER SPECTRUM (reciprocal space) ======= + ax2 = nexttile; + imagesc(kx, ky, log(1 + abs(fft_imgs{k}).^2)); + axis image; + set(gca, 'FontSize', 14, 'YDir', 'normal') + xlabel('k_x [\mum^{-1}]', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font); + ylabel('k_y [\mum^{-1}]', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font); + title('Power Spectrum - S(k_x,k_y)', 'Interpreter', 'tex', ... + 'FontSize', 16, 'FontWeight', 'bold', 'FontName', font); + 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 + drawPSOverlays(kx, ky, r_min, r_max) + + % ======= RADIAL DISTRIBUTION (S(k)) ======= + 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); + set(gca, 'FontSize', 14, 'YScale', 'log', 'XLim', [min(k_rho_vals), max(k_rho_vals)]); + xlabel('k_\rho [\mum^{-1}]', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font); + ylabel('Magnitude (a.u.)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font); + title('Radial Spectral Distribution - S(k_\rho)', 'Interpreter', 'tex', ... + '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) + % ======= ANGULAR DISTRIBUTION (S(θ)) ======= + nexttile; + plot(theta_vals/pi, S_theta, 'LineWidth', 2); + set(gca, 'FontSize', 14, 'YScale', 'log', 'YLim', [1E4, 1E7]); + xlabel('\theta/\pi [rad]', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font); + ylabel('Magnitude (a.u.)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font); + title('Angular Spectral Distribution - S(\theta)', 'Interpreter', 'tex', ... + 'FontSize', 16, 'FontWeight', 'bold', 'FontName', font); + grid on; % Enable major grid + ax = gca; + ax.MinorGridLineStyle = ':'; % Optional: make minor grid dotted + ax.MinorGridColor = [0.7 0.7 0.7]; % Optional: light gray minor grid color + ax.MinorGridAlpha = 0.5; % Optional: transparency for minor grid - % 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 + ax.XMinorGrid = 'on'; % Enable minor grid for x-axis + ax.YMinorGrid = 'on'; % Enable minor grid for y-axis (if desired) + + drawnow; + + if ~fft_size_known + fft_sz = size(fft_imgs{k}); + N_radial_bins_used = length(S_k_smoothed); + N_angular_bins_used = length(S_theta); + avg_ps_accum = zeros(fft_sz); + avg_S_k_accum = zeros(1, N_radial_bins_used); + avg_S_theta_accum = zeros(1, N_angular_bins_used); + fft_size_known = true; + end + + avg_ps_accum = avg_ps_accum + abs(fft_imgs{k}).^2; + avg_S_k_accum = avg_S_k_accum + S_k_smoothed; + avg_S_theta_accum = avg_S_theta_accum + S_theta; + + if ~skipMovieRender + % 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 + if ~skipSaveFigures + % Construct a filename for each image + fileNamePNG = fullfile(saveFolder, sprintf('fft_analysis_img_%03d.png', k)); + + % Save current figure as PNG with high resolution + print(gcf, fileNamePNG, '-dpng', '-r100'); % 300 dpi for high quality + end + if skipMovieRender & skipSaveFigures + pause(0.5); + end end -% Close the video file -close(videoFile); +if ~skipMovieRender + % Close the video file + close(videoFile); +end + +%% ===== Final Averages ===== +avg_ps = avg_ps_accum / N_shots; +avg_S_k = avg_S_k_accum / N_shots; +avg_S_theta = avg_S_theta_accum / N_shots; + +% Generate figure with 3 subplots +figure('Name', 'Average Spectral Analysis', 'Position', [400 200 1200 400]); +tavg = tiledlayout(1, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); + +% ==== 1. Average FFT Power Spectrum ==== +nexttile; +imagesc(kx, ky, log(1 + avg_ps)); +axis image; +set(gca, 'FontSize', 14, 'YDir', 'normal') +xlabel('k_x [\mum^{-1}]', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font); +ylabel('k_y [\mum^{-1}]', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font); +title('Average Power Spectrum', 'FontSize', 16, 'FontWeight', 'bold'); +colorbar; +colormap(Colormaps.inferno()); + +% ==== 2. Average Radial Spectral Distribution ==== +nexttile; +plot(k_rho_vals, avg_S_k, 'LineWidth', 2); +xlabel('k_\rho [\mum^{-1}]', 'Interpreter', 'tex', 'FontSize', 14); +ylabel('Magnitude (a.u.)', 'Interpreter', 'tex', 'FontSize', 14); +title('Average S(k_\rho)', 'FontSize', 16, 'FontWeight', 'bold'); +set(gca, 'FontSize', 14, 'YScale', 'log', 'XLim', [min(k_rho_vals), max(k_rho_vals)]); +grid on; + +% ==== 3. Average Angular Spectral Distribution ==== +nexttile; +plot(theta_vals/pi, avg_S_theta, 'LineWidth', 2); +xlabel('\theta/\pi [rad]', 'Interpreter', 'tex', 'FontSize', 14); +ylabel('Magnitude (a.u.)', 'Interpreter', 'tex', 'FontSize', 14); +title('Average S(\theta)', 'FontSize', 16, 'FontWeight', 'bold'); +set(gca, 'FontSize', 14, 'YScale', 'log'); +grid on; +ax = gca; +ax.XMinorGrid = 'on'; +ax.YMinorGrid = 'on'; %% Helper Functions function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization) @@ -382,55 +466,41 @@ function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMas 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) +function [k_rho_vals, S_radial] = computeRadialSpectralDistribution(IMGFFT, kx, ky, thetamin, thetamax, num_bins) + % IMGFFT : 2D FFT image (fftshifted and cropped) + % kx, ky : 1D physical wavenumber axes [μm⁻¹] matching FFT size + % thetamin : Minimum angle (in radians) + % thetamax : Maximum angle (in radians) + % num_bins : Number of radial 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; + [KX, KY] = meshgrid(kx, ky); + K_rho = sqrt(KX.^2 + KY.^2); + Theta = atan2(KY, KX); - % 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); + power_spectrum = abs(IMGFFT).^2; - % Radial bins + r_min = min(K_rho(angle_mask)); + r_max = max(K_rho(angle_mask)); 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); + radial_mask = (K_rho >= r_low) & (K_rho < 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) +function [theta_vals, S_theta] = computeAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma, windowSize) % Apply threshold to isolate strong peaks IMGFFT(IMGFFT < threshold) = 0; @@ -477,9 +547,6 @@ function [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IM 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) @@ -561,28 +628,235 @@ function ret = cropODImage(img, center, span) 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 +function imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime) +%CALCULATEODIMAGE Calculates the optical density (OD) image for absorption imaging. +% +% imageOD = calculateODImage(imageAtom, imageBackground, imageDark, mode, exposureTime) +% +% Inputs: +% imageAtom - Image with atoms +% imageBackground - Image without atoms +% imageDark - Image without light +% mode - 'LowIntensity' (default) or 'HighIntensity' +% exposureTime - Required only for 'HighIntensity' [in seconds] +% +% Output: +% imageOD - Computed OD image +% + arguments + imageAtom (:,:) {mustBeNumeric} + imageBackground (:,:) {mustBeNumeric} + imageDark (:,:) {mustBeNumeric} + mode char {mustBeMember(mode, {'LowIntensity', 'HighIntensity'})} = 'LowIntensity' + exposureTime double = NaN + end + + % Compute numerator and denominator numerator = imageBackground - imageDark; denominator = imageAtom - imageDark; + % Avoid division by zero numerator(numerator == 0) = 1; denominator(denominator == 0) = 1; - - ret = -log(double(abs(denominator ./ numerator))); - if numel(ret) == 1 - ret = ret(1); + % Calculate OD based on mode + switch mode + case 'LowIntensity' + imageOD = -log(abs(denominator ./ numerator)); + + case 'HighIntensity' + if isnan(exposureTime) + error('Exposure time must be provided for HighIntensity mode.'); + end + imageOD = abs(denominator ./ numerator); + imageOD = -log(imageOD) + (numerator - denominator) ./ (7000 * (exposureTime / 5e-6)); end + +end + +function drawODOverlays(x1, y1, x2, y2) + + % Parameters + tick_spacing = 10; % µm between ticks + tick_length = 2; % µm tick mark length + line_color = [0.5 0.5 0.5]; + tick_color = [0.5 0.5 0.5]; + font_size = 10; + + % Vector from start to end + dx = x2 - x1; + dy = y2 - y1; + L = sqrt(dx^2 + dy^2); + + % Unit direction vector along diagonal + ux = dx / L; + uy = dy / L; + + % Perpendicular unit vector for ticks + perp_ux = -uy; + perp_uy = ux; + + % Midpoint (center) + xc = (x1 + x2) / 2; + yc = (y1 + y2) / 2; + + % Number of positive and negative ticks + n_ticks = floor(L / (2 * tick_spacing)); + + % Draw main diagonal line + plot([x1 x2], [y1 y2], '--', 'Color', line_color, 'LineWidth', 1.2); + + for i = -n_ticks:n_ticks + d = i * tick_spacing; + xt = xc + d * ux; + yt = yc + d * uy; + + % Tick line endpoints + xt1 = xt - 0.5 * tick_length * perp_ux; + yt1 = yt - 0.5 * tick_length * perp_uy; + xt2 = xt + 0.5 * tick_length * perp_ux; + yt2 = yt + 0.5 * tick_length * perp_uy; + + % Draw tick + plot([xt1 xt2], [yt1 yt2], '--', 'Color', tick_color, 'LineWidth', 1); + + % Label: centered at tick, offset slightly along diagonal + if d ~= 0 + text(xt, yt, sprintf('%+d', d), ... + 'Color', tick_color, ... + 'FontSize', font_size, ... + 'HorizontalAlignment', 'center', ... + 'VerticalAlignment', 'bottom', ... + 'Rotation', atan2d(dy, dx)); + end + + end +end + +function drawPSOverlays(kx, ky, r_min, r_max) +% drawFFTOverlays - Draw overlays on existing FFT plot: +% - Radial lines every 30° +% - Annular highlight with white (upper half) and gray (lower half) circles between r_min and r_max +% - Horizontal white bands at ky=0 in annulus region +% - Scale ticks and labels every 1 μm⁻¹ along each radial line +% +% Inputs: +% kx, ky - reciprocal space vectors (μm⁻¹) +% r_min - inner annulus radius offset index (integer) +% r_max - outer annulus radius offset index (integer) +% +% Example: +% hold on; +% drawFFTOverlays(kx, ky, 10, 30); + + hold on + + % === Overlay Radial Lines + Scales === + [kx_grid, ky_grid] = meshgrid(kx, ky); + [~, kr_grid] = cart2pol(kx_grid, ky_grid); % kr_grid in μm⁻¹ + + max_kx = max(kx); + max_ky = max(ky); + + for angle = 0 : pi/6 : pi + x_line = [0, max_kx] * cos(angle); + y_line = [0, max_ky] * sin(angle); + + % Plot radial lines + plot(x_line, y_line, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2); + plot(x_line, -y_line, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2); + + % Draw scale ticks along positive radial line + drawTicksAlongLine(0, 0, x_line(2), y_line(2)); + + % Draw scale ticks along negative radial line (reflect y) + drawTicksAlongLine(0, 0, x_line(2), -y_line(2)); + end + + % === Overlay Annular Highlight: White (r_min to r_max), Gray elsewhere === + theta_full = linspace(0, 2*pi, 500); + + center_x = ceil(size(kr_grid, 2) / 2); + center_y = ceil(size(kr_grid, 1) / 2); + + k_min = kr_grid(center_y, center_x + r_min); + k_max = kr_grid(center_y, center_x + r_max); + + % Upper half: white dashed circles + x1_upper = k_min * cos(theta_full(theta_full <= pi)); + y1_upper = k_min * sin(theta_full(theta_full <= pi)); + x2_upper = k_max * cos(theta_full(theta_full <= pi)); + y2_upper = k_max * sin(theta_full(theta_full <= pi)); + plot(x1_upper, y1_upper, 'w--', 'LineWidth', 1.2); + plot(x2_upper, y2_upper, 'w--', 'LineWidth', 1.2); + + % Lower half: gray dashed circles + x1_lower = k_min * cos(theta_full(theta_full > pi)); + y1_lower = k_min * sin(theta_full(theta_full > pi)); + x2_lower = k_max * cos(theta_full(theta_full > pi)); + y2_lower = k_max * sin(theta_full(theta_full > pi)); + plot(x1_lower, y1_lower, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.0); + plot(x2_lower, y2_lower, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.0); + + % === Highlight horizontal band across k_y = 0 === + x_vals = kx; + xW1 = x_vals((x_vals >= -k_max) & (x_vals < -k_min)); + xW2 = x_vals((x_vals > k_min) & (x_vals <= k_max)); + + plot(xW1, zeros(size(xW1)), 'w--', 'LineWidth', 1.2); + plot(xW2, zeros(size(xW2)), 'w--', 'LineWidth', 1.2); + + hold off + + + % --- Nested helper function to draw ticks along a radial line --- + function drawTicksAlongLine(x_start, y_start, x_end, y_end) + % Tick parameters + tick_spacing = 1; % spacing between ticks in μm⁻¹ + tick_length = 0.05 * sqrt((x_end - x_start)^2 + (y_end - y_start)^2); % relative tick length + line_color = [0.5 0.5 0.5]; + tick_color = [0.5 0.5 0.5]; + font_size = 8; + + % Vector along the line + dx = x_end - x_start; + dy = y_end - y_start; + L = sqrt(dx^2 + dy^2); + ux = dx / L; + uy = dy / L; + + % Perpendicular vector for ticks + perp_ux = -uy; + perp_uy = ux; + + % Number of ticks (from 0 up to max length) + n_ticks = floor(L / tick_spacing); + + for i = 1:n_ticks + % Position of tick along the line + xt = x_start + i * tick_spacing * ux; + yt = y_start + i * tick_spacing * uy; + + % Tick endpoints + xt1 = xt - 0.5 * tick_length * perp_ux; + yt1 = yt - 0.5 * tick_length * perp_uy; + xt2 = xt + 0.5 * tick_length * perp_ux; + yt2 = yt + 0.5 * tick_length * perp_uy; + + % Draw tick + plot([xt1 xt2], [yt1 yt2], '-', 'Color', tick_color, 'LineWidth', 1); + + % Label with distance (integer) + text(xt, yt, sprintf('%d', i), ... + 'Color', tick_color, ... + 'FontSize', font_size, ... + 'HorizontalAlignment', 'center', ... + 'VerticalAlignment', 'bottom', ... + 'Rotation', atan2d(dy, dx)); + end + end + end function [optrefimages] = removefringesInImage(absimages, refimages, bgmask) @@ -654,4 +928,4 @@ function [optrefimages] = removefringesInImage(absimages, refimages, bgmask) % Compute optimised reference image optrefimages(:,:,j)=reshape(R*c,[ydim xdim]); end -end +end \ No newline at end of file