MAJOR update - many additions to plotting
This commit is contained in:
parent
cd48c0ab9b
commit
cc72fb2ccc
@ -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
|
||||
t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact');
|
||||
|
||||
% 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
|
||||
% ======= 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
|
||||
imagesc(x, y, IMG)
|
||||
|
||||
% 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
|
||||
|
||||
% 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;
|
||||
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
|
||||
ylabel(hcb, 'Optical Density', 'Rotation', -90, 'FontSize', 14, 'FontName', font);
|
||||
|
||||
% Plot the smoothed radial distribution
|
||||
nexttile
|
||||
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());
|
||||
|
||||
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;
|
||||
|
||||
% ======= 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
|
||||
|
||||
% 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
|
||||
ax.XMinorGrid = 'on'; % Enable minor grid for x-axis
|
||||
ax.YMinorGrid = 'on'; % Enable minor grid for y-axis (if desired)
|
||||
|
||||
drawnow
|
||||
pause(0.5)
|
||||
drawnow;
|
||||
|
||||
% 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
|
||||
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)));
|
||||
% Calculate OD based on mode
|
||||
switch mode
|
||||
case 'LowIntensity'
|
||||
imageOD = -log(abs(denominator ./ numerator));
|
||||
|
||||
if numel(ret) == 1
|
||||
ret = ret(1);
|
||||
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)
|
||||
|
||||
Loading…
Reference in New Issue
Block a user