197 lines
7.3 KiB
Matlab
197 lines
7.3 KiB
Matlab
function [UCD] = extractAverageUnitCellDensity(folder_path, run_index, radius, minPeakHeight, SuppressPlotFlag)
|
|
|
|
set(0,'defaulttextInterpreter','latex')
|
|
set(groot, 'defaultAxesTickLabelInterpreter','latex');
|
|
set(groot, 'defaultLegendInterpreter','latex');
|
|
|
|
% Load data
|
|
filePath = sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'), run_index);
|
|
|
|
try
|
|
Data = load(filePath, 'psi', 'Params', 'Transf', 'Observ');
|
|
catch ME
|
|
warning('Failed to load file: %s\n%s', filePath, ME.message);
|
|
UCD = NaN;
|
|
return;
|
|
end
|
|
|
|
Params = Data.Params;
|
|
Transf = Data.Transf;
|
|
psi = gather(Data.psi);
|
|
|
|
% Axes in micrometers
|
|
x = Transf.x * Params.l0 * 1e6;
|
|
y = Transf.y * Params.l0 * 1e6;
|
|
z = Transf.z * Params.l0 * 1e6;
|
|
dz = z(2) - z(1);
|
|
|
|
% Compute integrated density
|
|
n = abs(psi).^2;
|
|
nxy = squeeze(trapz(n * dz, 3));
|
|
state = nxy';
|
|
|
|
% Fourier spectrum for orientation detection
|
|
F = fftshift(abs(fft2(state - mean(state(:)))));
|
|
[nx, ny] = size(state);
|
|
kx = linspace(-pi / (x(2) - x(1)), pi / (x(2) - x(1)), nx);
|
|
ky = linspace(-pi / (y(2) - y(1)), pi / (y(2) - y(1)), ny);
|
|
[KX, KY] = meshgrid(kx, ky);
|
|
theta = atan2(KY, KX);
|
|
|
|
% Remove center/DC
|
|
R = sqrt(KX.^2 + KY.^2);
|
|
F(R < 0.1 * max(kx(:))) = 0;
|
|
|
|
% Angular histogram of power spectrum
|
|
nbins = 180;
|
|
angles = linspace(-pi, pi, nbins);
|
|
powerAngular = zeros(1, nbins);
|
|
for k = 1:nbins-1
|
|
mask = theta >= angles(k) & theta < angles(k+1);
|
|
powerAngular(k) = sum(F(mask), 'all');
|
|
end
|
|
powerAngular(end) = powerAngular(1); % wrap around
|
|
|
|
% Anisotropy measure
|
|
peakRatio = max(powerAngular) / mean(powerAngular);
|
|
|
|
% Threshold to distinguish stripe vs droplet
|
|
isStripe = peakRatio > 20;
|
|
|
|
if ~isStripe
|
|
% === DROPLET MODE: Voronoi cell ===
|
|
peaks = extractPeaks(state, radius, minPeakHeight);
|
|
[row, col] = find(peaks);
|
|
peakHeights = state(peaks);
|
|
[~, maxIdx] = max(peakHeights);
|
|
MaxPeakLocation = [row(maxIdx), col(maxIdx)];
|
|
points = [col, row];
|
|
|
|
[V, C] = voronoin(points);
|
|
Vcell_indices = C{maxIdx};
|
|
|
|
if all(Vcell_indices > 0) && all(Vcell_indices <= size(V,1)) && ~any(isinf(V(Vcell_indices, 1)))
|
|
vx = interp1(1:numel(x), x, V(Vcell_indices,1), 'linear', 'extrap');
|
|
vy = interp1(1:numel(y), y, V(Vcell_indices,2), 'linear', 'extrap');
|
|
vx(end+1) = vx(1);
|
|
vy(end+1) = vy(1);
|
|
AreaOfVoronoiCell = polyarea(vx, vy);
|
|
|
|
[X, Y] = meshgrid(x, y);
|
|
inCellMask = inpolygon(X, Y, vx, vy);
|
|
NumberOfParticles = sum(state(inCellMask));
|
|
UCD = NumberOfParticles / AreaOfVoronoiCell;
|
|
else
|
|
warning('Voronoi cell for max peak is invalid.');
|
|
UCD = NaN;
|
|
end
|
|
|
|
else
|
|
% === STRIPE MODE: Use rectangular unit cell aligned with stripe direction ===
|
|
[~, idxMax] = max(F(:));
|
|
stripeAngle = theta(idxMax); % radians
|
|
|
|
% Rotate image to align stripes horizontally
|
|
rotatedState = imrotate(state, -rad2deg(stripeAngle), 'bilinear', 'crop');
|
|
|
|
% Estimate vertical stripe spacing (stripe-normal direction)
|
|
stripeProfileY = mean(rotatedState, 2) - mean(rotatedState(:));
|
|
fftY = abs(fft(stripeProfileY));
|
|
fftY(1:3) = 0;
|
|
[~, fyIdx] = max(fftY(1:floor(end/2)));
|
|
spacingY = round(numel(stripeProfileY) / fyIdx);
|
|
|
|
% Estimate horizontal spacing (along-stripe periodicity)
|
|
stripeProfileX = mean(rotatedState, 1) - mean(rotatedState(:));
|
|
fftX = abs(fft(stripeProfileX));
|
|
fftX(1:3) = 0;
|
|
[~, fxIdx] = max(fftX(1:floor(end/2)));
|
|
spacingX = round(numel(stripeProfileX) / fxIdx);
|
|
|
|
% Find stripe center (vertical)
|
|
[~, centerY] = max(stripeProfileY);
|
|
rowRange = max(1, centerY - round(0.5 * spacingY)) : ...
|
|
min(size(rotatedState,1), centerY + round(0.5 * spacingY));
|
|
|
|
% Find repeat center along stripe direction (horizontal)
|
|
[~, centerX] = max(stripeProfileX);
|
|
colRange = max(1, centerX - round(0.5 * spacingX)) : ...
|
|
min(size(rotatedState,2), centerX + round(0.5 * spacingX));
|
|
|
|
% Extract unit cell region
|
|
unitCellRegion = rotatedState(rowRange, colRange);
|
|
NumberOfParticles = sum(unitCellRegion(:));
|
|
AreaOfUnitCell = numel(unitCellRegion) * abs(x(2)-x(1)) * abs(y(2)-y(1));
|
|
UCD = NumberOfParticles / AreaOfUnitCell;
|
|
end
|
|
|
|
|
|
% Optional plot
|
|
if ~SuppressPlotFlag
|
|
figure(1);
|
|
clf
|
|
set(gcf,'Position', [100, 100, 900, 900])
|
|
if ~isStripe
|
|
plotxy = pcolor(x,y,state);
|
|
set(plotxy, 'EdgeColor', 'none');
|
|
hold on;
|
|
plot(vx, vy, 'w-', 'LineWidth', 2);
|
|
cbar1 = colorbar;
|
|
cbar1.Label.Interpreter = 'latex';
|
|
% cbar1.Ticks = []; % Disable the ticks
|
|
colormap(gca, Helper.Colormaps.plasma())
|
|
xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16)
|
|
ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16)
|
|
title(['UCD = ' num2str(UCD) ' \mum^{-2}'], ...
|
|
'Interpreter', 'tex', 'FontSize', 16)
|
|
set(gca, 'FontSize', 16); % For tick labels only
|
|
else
|
|
imagesc(rotatedState);
|
|
axis image;
|
|
cbar1 = colorbar;
|
|
cbar1.Label.Interpreter = 'latex';
|
|
colormap(gca, Helper.Colormaps.plasma())
|
|
xlabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16)
|
|
ylabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 16)
|
|
|
|
% Title with UCD
|
|
title(['UCD = ' num2str(UCD, '%.3f') ' $\mu$m$^{-2}$'], ...
|
|
'Interpreter', 'latex', 'FontSize', 16)
|
|
|
|
% Highlight unit cell region
|
|
rectX = colRange(1) - 0.5;
|
|
rectY = rowRange(1) - 0.5;
|
|
rectW = length(colRange);
|
|
rectH = length(rowRange);
|
|
rectangle('Position', [rectX, rectY, rectW, rectH], ...
|
|
'EdgeColor', 'w', 'LineWidth', 2);
|
|
title(sprintf('UCD = %.4f $\\mu$m$^{-2}$', UCD), ...
|
|
'Interpreter', 'latex', 'FontSize', 16);
|
|
|
|
|
|
end
|
|
end
|
|
end
|
|
|
|
function [peaks_mask] = extractPeaks(image, radius, minPeakHeight)
|
|
peaks = imregionalmax(image);
|
|
peaks(image < minPeakHeight) = 0;
|
|
|
|
[row, col] = find(peaks); % row and col are the coordinates of the detected peaks
|
|
unique_peaks = true(size(row)); % Assume all peaks are unique initially
|
|
|
|
% Check distances between each peak and remove duplicates
|
|
for i = 1:length(row)
|
|
if unique_peaks(i) % If this peak hasn't been excluded
|
|
% Find the distance from this peak to all other peaks
|
|
dist = sqrt((row(i) - row).^2 + (col(i) - col).^2); % Euclidean distance
|
|
|
|
% Exclude any peak within the radius
|
|
unique_peaks(dist < radius & dist > 0) = false; % Exclude other peaks in radius
|
|
end
|
|
end
|
|
|
|
% Save the peaks
|
|
peaks_mask = false(size(image));
|
|
peaks_mask(sub2ind(size(image), row(unique_peaks), col(unique_peaks))) = true;
|
|
end |