99 lines
3.3 KiB
Matlab
99 lines
3.3 KiB
Matlab
function extractHexaticOrder(I)
|
|
% Input:
|
|
% I = 2D image matrix (grayscale, double)
|
|
% Output:
|
|
% Plots original image + hexatic order + histogram
|
|
|
|
%% 1) Detect dots (local maxima)
|
|
I_smooth = imgaussfilt(I, 2);
|
|
bw = imregionalmax(I_smooth);
|
|
[y, x] = find(bw);
|
|
|
|
if length(x) < 6
|
|
error('Too few detected dots (%d) for meaningful hexatic order.', length(x));
|
|
end
|
|
|
|
%% 2) Compute hexatic order parameter psi6
|
|
|
|
% Delaunay triangulation to find neighbors
|
|
tri = delaunay(x,y);
|
|
N = length(x);
|
|
neighbors = cell(N,1);
|
|
for i = 1:size(tri,1)
|
|
neighbors{tri(i,1)} = unique([neighbors{tri(i,1)} tri(i,2) tri(i,3)]);
|
|
neighbors{tri(i,2)} = unique([neighbors{tri(i,2)} tri(i,1) tri(i,3)]);
|
|
neighbors{tri(i,3)} = unique([neighbors{tri(i,3)} tri(i,1) tri(i,2)]);
|
|
end
|
|
|
|
psi6_vals = zeros(N,1);
|
|
for j = 1:N
|
|
nbrs = neighbors{j};
|
|
nbrs(nbrs == j) = [];
|
|
if isempty(nbrs)
|
|
psi6_vals(j) = 0;
|
|
continue;
|
|
end
|
|
angles = atan2(y(nbrs)-y(j), x(nbrs)-x(j));
|
|
psi6_vals(j) = abs(mean(exp(1i*6*angles)));
|
|
end
|
|
psi6_global = mean(psi6_vals);
|
|
|
|
%% 3) Plotting
|
|
|
|
figure('Name','Hexatic Order Analysis','NumberTitle','off', 'Units','normalized', 'Position',[0.2 0.2 0.6 0.7]);
|
|
t = tiledlayout(2,2, 'Padding','none', 'TileSpacing','compact');
|
|
|
|
% Top left: Original image with detected dots (square)
|
|
ax1 = nexttile(1);
|
|
imshow(I, []);
|
|
hold on;
|
|
scatter(x, y, 30, 'w', 'filled');
|
|
hold off;
|
|
colormap(ax1, Helper.Colormaps.plasma()); % Example colormap for original image plot
|
|
cbar1 = colorbar(ax1);
|
|
cbar1.Label.Interpreter = 'latex';
|
|
xlabel(ax1, '$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14);
|
|
ylabel(ax1, '$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14);
|
|
title(ax1, '$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14);
|
|
axis(ax1, 'square');
|
|
|
|
% Top right: Histogram of local |psi6| (square)
|
|
ax2 = nexttile(2);
|
|
histogram(psi6_vals, 20, 'FaceColor', 'b');
|
|
xlabel(ax2, '$|\psi_6|$', 'Interpreter', 'latex', 'FontSize', 14);
|
|
ylabel(ax2, 'Count');
|
|
title(ax2, 'Histogram of Local Hexatic Order');
|
|
xlim(ax2, [0 1]);
|
|
grid(ax2, 'on');
|
|
axis(ax2, 'square');
|
|
|
|
% Bottom: Hexatic order magnitude on dots with bond orientation arrows (span 2 tiles)
|
|
ax3 = nexttile([1 2]);
|
|
scatter_handle = scatter(x, y, 50, psi6_vals, 'filled');
|
|
colormap(ax3, 'jet'); % Different colormap for hexatic order magnitude
|
|
cbar3 = colorbar(ax3);
|
|
cbar3.Label.String = '|$\psi_6$|';
|
|
cbar3.Label.Interpreter = 'latex';
|
|
caxis(ax3, [0 1]);
|
|
axis(ax3, 'equal');
|
|
axis(ax3, 'tight');
|
|
title(ax3, sprintf('Local Hexatic Order |\\psi_6|, Global = %.3f', psi6_global));
|
|
hold(ax3, 'on');
|
|
|
|
for j = 1:N
|
|
nbrs = neighbors{j};
|
|
nbrs(nbrs == j) = [];
|
|
if isempty(nbrs)
|
|
continue;
|
|
end
|
|
angles = atan2(y(nbrs)-y(j), x(nbrs)-x(j));
|
|
psi6_complex = mean(exp(1i*6*angles));
|
|
local_angle = angle(psi6_complex)/6;
|
|
arrow_length = 5;
|
|
quiver(ax3, x(j), y(j), arrow_length*cos(local_angle), arrow_length*sin(local_angle), ...
|
|
'k', 'MaxHeadSize', 2, 'AutoScale', 'off');
|
|
end
|
|
hold(ax3, 'off');
|
|
|
|
end
|