Calculations/Data-Analyzer/conductSpectralAnalysis.m

1360 lines
47 KiB
Matlab
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

%% ===== D-S Settings =====
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 = "//DyLabNAS/Data/TwoDGas/2025/06/23/";
run = '0300';
folderPath = strcat(folderPath, run);
cam = 5;
angle = 0;
center = [1410, 2030];
span = [200, 200];
fraction = [0.1, 0.1];
pixel_size = 5.86e-6; % in meters
magnification = 23.94;
removeFringes = false;
ImagingMode = 'HighIntensity';
PulseDuration = 5e-6; % in s
% 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
% Plotting and saving
scan_parameter = 'ps_rot_mag_fin_pol_angle';
% scan_parameter = 'rot_mag_field';
savefileName = 'DropletsToStripes';
font = 'Bahnschrift';
if strcmp(savefileName, 'DropletsToStripes')
scan_groups = 0:5:45;
titleString = 'Droplets to Stripes';
elseif strcmp(savefileName, 'StripesToDroplets')
scan_groups = 45:-5:0;
titleString = 'Stripes to Droplets';
end
% Flags
skipNormalization = false;
skipUnshuffling = true;
skipPreprocessing = true;
skipMasking = true;
skipIntensityThresholding = true;
skipBinarization = true;
skipMovieRender = true;
skipSaveFigures = false;
skipSaveOD = false;
%% ===== S-D Settings =====
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 = "//DyLabNAS/Data/TwoDGas/2025/06/24/";
run = '0001';
folderPath = strcat(folderPath, run);
cam = 5;
angle = 0;
center = [1410, 2030];
span = [200, 200];
fraction = [0.1, 0.1];
pixel_size = 5.86e-6; % in meters
magnification = 23.94;
removeFringes = false;
ImagingMode = 'HighIntensity';
PulseDuration = 5e-6; % in s
% 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
% Plotting and saving
scan_parameter = 'ps_rot_mag_fin_pol_angle';
% scan_parameter = 'rot_mag_field';
savefileName = 'StripesToDroplets';
font = 'Bahnschrift';
if strcmp(savefileName, 'DropletsToStripes')
scan_groups = 0:5:45;
titleString = 'Droplets to Stripes';
elseif strcmp(savefileName, 'StripesToDroplets')
scan_groups = 45:-5:0;
titleString = 'Stripes to Droplets';
end
% Flags
skipNormalization = true;
skipUnshuffling = false;
skipPreprocessing = true;
skipMasking = true;
skipIntensityThresholding = true;
skipBinarization = true;
skipMovieRender = true;
skipSaveFigures = false;
skipSaveOD = false;
%% ===== Load and 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, ImagingMode, PulseDuration), 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);
radial_spectral_contrast = zeros(1, nimgs);
angular_spectral_weight = zeros(1, nimgs);
N_shots = length(od_imgs);
if ~skipMovieRender
% Create VideoWriter object for movie
videoFile = VideoWriter([savefileName '.mp4'], 'MPEG-4');
videoFile.Quality = 100; % Set quality to maximum (0100)
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
ps_list = cell(1, N_shots); % 2D power spectrum
s_k_list = cell(1, N_shots); % Radial spectrum
s_theta_list = cell(1, N_shots); % Angular spectrum
for k = 1:N_shots
IMG = od_imgs{k};
[IMGFFT, IMGPR] = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization);
% Size of original image (in pixels)
[Ny, Nx] = size(IMG);
% 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);
% 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
radial_spectral_contrast(k) = computeRadialSpectralContrast(fft_imgs{k}, r_min, r_max, Angular_Threshold);
S_theta_norm = S_theta / max(S_theta); % Normalize to 1
angular_spectral_weight(k) = trapz(theta_vals, S_theta_norm);
ps_list{k} = abs(fft_imgs{k}).^2; % store the power spectrum
s_k_list{k} = S_k_smoothed; % store smoothed radial spectrum
s_theta_list{k} = S_theta; % store angular spectrum
figure(1);
clf
set(gcf,'Position',[500 100 1000 800])
t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact');
% ======= OD IMAGE (real space) =======
ax1 = nexttile;
imagesc(x, y, IMG)
hold on;
% 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;
axis equal tight;
set(gca, 'FontSize', 14, 'YDir', 'normal')
colormap(ax1, Colormaps.inferno());
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);
if strcmp(scan_parameter, 'ps_rot_mag_fin_pol_angle')
text(0.975, 0.975, [num2str(scan_parameter_values(k), '%.1f^\\circ')], ...
'Color', 'white', 'FontWeight', 'bold', 'FontSize', 14, ...
'Interpreter', 'tex', 'Units', 'normalized', ...
'HorizontalAlignment', 'right', 'VerticalAlignment', 'top');
else
text(0.975, 0.975, [num2str(scan_parameter_values(k), '%.2f'), ' G'], ...
'Color', 'white', 'FontWeight', 'bold', 'FontSize', 14, ...
'Interpreter', 'tex', 'Units', 'normalized', ...
'HorizontalAlignment', 'right', 'VerticalAlignment', 'top');
end
% ======= FFT POWER SPECTRUM (reciprocal space) =======
ax2 = nexttile;
imagesc(kx, ky, log(1 + ps_list{k}));
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.coolwarm());
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, '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;
if ~skipNormalization
plot(theta_vals/pi, S_theta_norm, 'LineWidth', 2);
set(gca, 'FontSize', 14, 'YLim', [0, 1]);
else
plot(theta_vals/pi, S_theta, 'LineWidth', 2);
set(gca, 'FontSize', 14, 'YScale', 'log', 'YLim', [1E4, 1E7]);
end
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
ax.XMinorGrid = 'on'; % Enable minor grid for x-axis
ax.YMinorGrid = 'on'; % Enable minor grid for y-axis (if desired)
drawnow;
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 ~skipSaveOD
odDataStruct = struct();
odDataStruct.IMG = IMG;
odDataStruct.x = x;
odDataStruct.y = y;
odDataStruct.scan_parameter_value = scan_parameter_values(k);
save(fullfile(saveFolder, sprintf('od_image_%03d.mat', k)), '-struct', 'odDataStruct');
end
if skipMovieRender & skipSaveFigures
pause(0.5);
end
end
if ~skipMovieRender
% Close the video file
close(videoFile);
end
%% Track across the transition
% 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_sc = zeros(size(unique_scan_parameter_values));
stderr_sc = zeros(size(unique_scan_parameter_values));
% Loop through each unique theta and compute mean and standard error
for i = 1:length(unique_scan_parameter_values)
group_vals = radial_spectral_contrast(idx == i);
mean_sc(i) = mean(group_vals);
stderr_sc(i) = std(group_vals) / sqrt(length(group_vals)); % standard error = std / sqrt(N)
end
figure(2);
set(gcf,'Position',[100 100 950 750])
errorbar(unique_scan_parameter_values, mean_sc, stderr_sc, 'o--', ...
'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5);
set(gca, 'FontSize', 14); % For tick labels only
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
hYLabel = ylabel('Radial Spectral Contrast', 'Interpreter', 'tex');
hTitle = title(titleString, '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
grid on
% Preallocate arrays
mean_sw = zeros(size(unique_scan_parameter_values));
stderr_sw = zeros(size(unique_scan_parameter_values));
% Loop through each unique theta and compute mean and standard error
for i = 1:length(unique_scan_parameter_values)
group_vals = angular_spectral_weight(idx == i);
mean_sw(i) = mean(group_vals);
stderr_sw(i) = std(group_vals) / sqrt(length(group_vals)); % standard error = std / sqrt(N)
end
figure(3);
set(gcf,'Position',[100 100 950 750])
errorbar(unique_scan_parameter_values, mean_sw, stderr_sw, 'o--', ...
'LineWidth', 1.5, 'MarkerSize', 6, 'CapSize', 5);
set(gca, 'FontSize', 14); % For tick labels only
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
hYLabel = ylabel('Angular Spectral Weight', 'Interpreter', 'tex');
hTitle = title(titleString, '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
grid on;
%% Plot Averages
% Group by scan parameter values
[unique_scan_parameter_values, ~, idx] = unique(scan_parameter_values);
N_params = numel(unique_scan_parameter_values);
if ~skipSaveFigures
% Define folder for saving images
saveFolder = [savefileName '_SavedFigures'];
if ~exist(saveFolder, 'dir')
mkdir(saveFolder);
end
end
% Loop over each unique parameter value
for p = 1:N_params
current_param = unique_scan_parameter_values(p);
indices = find(idx == p); % Indices of shots for this parameter
N_shots = numel(indices);
% Initialize accumulators
avg_ps = 0;
avg_S_k = 0;
avg_S_theta = 0;
% Accumulate values
for j = 1:N_shots
avg_ps = avg_ps + ps_list{indices(j)};
avg_S_k = avg_S_k + s_k_list{indices(j)};
avg_S_theta = avg_S_theta + s_theta_list{indices(j)};
end
% Average over repetitions
avg_ps = avg_ps / N_shots;
avg_S_k = avg_S_k / N_shots;
avg_S_theta = avg_S_theta / N_shots;
% ==== Plot ====
figure(3);
set(gcf,'Position',[400 200 1200 400])
tavg = tiledlayout(1, 3, 'TileSpacing', 'compact', 'Padding', 'compact');
% 1. 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.coolwarm());
if strcmp(scan_parameter, 'ps_rot_mag_fin_pol_angle')
text(0.975, 0.975, [num2str(current_param, '%.1f^\\circ')], ...
'Color', 'white', 'FontWeight', 'bold', 'FontSize', 14, ...
'Interpreter', 'tex', 'Units', 'normalized', ...
'HorizontalAlignment', 'right', 'VerticalAlignment', 'top');
else
text(0.975, 0.975, [num2str(current_param, '%.2f'), ' G'], ...
'Color', 'white', 'FontWeight', 'bold', 'FontSize', 14, ...
'Interpreter', 'tex', 'Units', 'normalized', ...
'HorizontalAlignment', 'right', 'VerticalAlignment', 'top');
end
% 2. Radial Spectrum
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. Angular Spectrum
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', ...
'YLim', [1E4, 1E7]);
grid on;
ax = gca;
ax.XMinorGrid = 'on';
ax.YMinorGrid = 'on';
drawnow;
% ==== Save Figure ====
if ~skipSaveFigures
% Create a filename for each averaged plot
fileNamePNG = fullfile(saveFolder, sprintf('fft_avg_analysis_param_%03d.png', p));
% Save current figure as PNG with high resolution
print(gcf, fileNamePNG, '-dpng', '-r300'); % 300 dpi for high quality
else
pause(0.5)
end
end
%% ========= Replot OD images ==========
% Settings
filePattern = fullfile(saveFolder, 'od_image_*.mat');
files = dir(filePattern);
colormapName = 'inferno';
showText = true;
showOverlay = true;
font = 'Bahnschrift';
% Load and organize all OD images by parameter and repetition
nFiles = length(files);
if nFiles == 0
error('No .mat OD image files found in folder: %s', saveFolder);
end
% Load all data and extract parameter values
scan_values = zeros(1, nFiles);
allData = cell(1, nFiles);
for k = 1:nFiles
S = load(fullfile(files(k).folder, files(k).name));
scan_values(k) = S.scan_parameter_value;
allData{k} = S;
end
% Get unique parameter values
unique_params = unique(scan_values);
nParams = numel(unique_params);
% Group images: paramData{i} = [rep1, rep2, ...]
if strcmp(savefileName, 'StripesToDroplets')
unique_params = fliplr(unique_params);
end
paramData = cell(1, nParams);
for i = 1:nParams
idxs = find(scan_values == unique_params(i));
paramData{i} = allData(idxs);
end
% Get number of repetitions (assumes all same)
nReps = max(cellfun(@numel, paramData));
% Initialize figure with one row, nParams columns
figure(100); clf;
% Set number of columns (e.g., 4 or auto-compute from nParams)
nCols = min(4, nParams);
nRows = ceil(nParams / nCols);
% Create tiled layout with multiple rows
t = tiledlayout(nRows, nCols, 'TileSpacing', 'compact', 'Padding', 'compact');
% Adjust figure size accordingly
set(gcf, 'Position', [100 100 300*nCols 300*nRows]);
% Pre-create image handles
axesArray = gobjects(1, nParams);
imgArray = gobjects(1, nParams);
textArray = gobjects(1, nParams);
for i = 1:nParams
S = paramData{i}{1}; % First repetition to initialize
ax = nexttile(i);
axesArray(i) = ax;
imgArray(i) = imagesc(S.x, S.y, S.IMG);
axis equal tight;
set(ax, 'YDir', 'normal');
colormap(ax, Colormaps.(colormapName)());
colorbar;
xlabel('x [\mum]', 'Interpreter', 'tex', 'FontSize', 12, 'FontName', font);
ylabel('y [\mum]', 'Interpreter', 'tex', 'FontSize', 12, 'FontName', font);
if showOverlay
hold on;
drawODOverlays(S.x(1), S.y(1), S.x(end), S.y(end));
drawODOverlays(S.x(end), S.y(1), S.x(1), S.y(end));
hold off;
end
% Add initial label
if showText
if strcmp(scan_parameter, 'ps_rot_mag_fin_pol_angle')
labelStr = sprintf('%.1f^\\circ', S.scan_parameter_value);
else
labelStr = sprintf('%.2f G', S.scan_parameter_value);
end
textArray(i) = text(ax, 0.975, 0.975, labelStr, ...
'Color', 'white', 'FontWeight', 'bold', ...
'FontSize', 12, 'Interpreter', 'tex', ...
'Units', 'normalized', ...
'HorizontalAlignment', 'right', ...
'VerticalAlignment', 'top');
end
end
% 🔁 Loop over repetitions
for rep = 1:nReps
for i = 1:nParams
repsForParam = paramData{i};
if rep <= numel(repsForParam)
S = repsForParam{rep};
imgArray(i).CData = S.IMG;
% Update text if needed (optional)
if showText
if strcmp(scan_parameter, 'ps_rot_mag_fin_pol_angle')
labelStr = sprintf('%.1f^\\circ', S.scan_parameter_value);
else
labelStr = sprintf('%.2f G', S.scan_parameter_value);
end
textArray(i).String = labelStr;
end
end
end
drawnow; % Update figure
% Optional: pause or save frame
pause(0.2);
end
%% ========= Replot OD images in chunks by parameter ==========
% Settings
filePattern = fullfile(saveFolder, 'od_image_*.mat');
files = dir(filePattern);
colormapName = 'inferno';
showText = true;
showOverlay = true;
font = 'Bahnschrift';
paramStep = 2; % Show every paramStep-th parameter
pauseTime = 0.2; % Seconds between repetitions
% Load and organize all OD images
nFiles = numel(files);
scan_values = zeros(1, nFiles);
allData = cell(1, nFiles);
for k = 1:nFiles
S = load(fullfile(files(k).folder, files(k).name));
scan_values(k) = S.scan_parameter_value;
allData{k} = S;
end
% Sort and group by unique parameter values
[unique_params, ~, ic] = unique(scan_values);
nParams = numel(unique_params);
paramGroups = cell(1, nParams);
for i = 1:nParams
paramGroups{i} = allData(ic == i);
end
if strcmp(savefileName, 'StripesToDroplets')
unique_params = fliplr(unique_params);
paramGroups = fliplr(paramGroups);
end
% Select a subset of parameters
selectedIdx = 1:paramStep:nParams;
nDisplayParams = numel(selectedIdx);
selectedGroups = paramGroups(selectedIdx);
% Get max number of repetitions
nReps = max(cellfun(@numel, selectedGroups));
% Initialize figure
figure(101); clf;
tiledlayout(1, nDisplayParams, 'TileSpacing', 'compact', 'Padding', 'compact');
set(gcf, 'Position', [100 100 300*nDisplayParams 300]);
imgArray = gobjects(1, nDisplayParams);
textArray = gobjects(1, nDisplayParams);
% Initial plot (repetition 1)
for j = 1:nDisplayParams
ax = nexttile;
group = selectedGroups{j};
S = group{1};
imgArray(j) = imagesc(S.x, S.y, S.IMG);
axis equal tight;
set(ax, 'YDir', 'normal');
colormap(ax, Colormaps.(colormapName)());
colorbar;
xlabel('x [\mum]', 'Interpreter', 'tex', 'FontSize', 12, 'FontName', font);
ylabel('y [\mum]', 'Interpreter', 'tex', 'FontSize', 12, 'FontName', font);
if showOverlay
hold on;
drawODOverlays(S.x(1), S.y(1), S.x(end), S.y(end));
drawODOverlays(S.x(end), S.y(1), S.x(1), S.y(end));
hold off;
end
if showText
if strcmp(scan_parameter, 'ps_rot_mag_fin_pol_angle')
labelStr = sprintf('%.1f^\\circ', S.scan_parameter_value);
else
labelStr = sprintf('%.2f G', S.scan_parameter_value);
end
textArray(j) = text(0.975, 0.975, labelStr, ...
'Color', 'white', 'FontWeight', 'bold', ...
'FontSize', 12, 'Interpreter', 'tex', ...
'Units', 'normalized', ...
'HorizontalAlignment', 'right', ...
'VerticalAlignment', 'top');
end
end
% Loop through repetitions
for rep = 1:nReps
for j = 1:nDisplayParams
group = selectedGroups{j};
if rep <= numel(group)
S = group{rep};
imgArray(j).CData = S.IMG;
if showText
if strcmp(scan_parameter, 'ps_rot_mag_fin_pol_angle')
textArray(j).String = sprintf('%.1f^\\circ', S.scan_parameter_value);
else
textArray(j).String = sprintf('%.2f G', S.scan_parameter_value);
end
end
end
end
drawnow;
pause(pauseTime);
end
%% 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, 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
[KX, KY] = meshgrid(kx, ky);
K_rho = sqrt(KX.^2 + KY.^2);
Theta = atan2(KY, KX);
if thetamin < thetamax
angle_mask = (Theta >= thetamin) & (Theta <= thetamax);
else
angle_mask = (Theta >= thetamin) | (Theta <= thetamax);
end
power_spectrum = abs(IMGFFT).^2;
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);
for i = 1:num_bins
r_low = r_edges(i);
r_high = r_edges(i + 1);
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] = computeAngularSpectralDistribution(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
end
function contrast = computeRadialSpectralContrast(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 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;
% 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, 'k--', 'LineWidth', 1.2);
plot(x2_upper, y2_upper, 'k--', '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)), 'k--', 'LineWidth', 1.2);
plot(xW2, zeros(size(xW2)), 'k--', '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)
% 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