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