function [psi_global, order_ratio, Npoints] = computeBondOrderParameters(I)
% Analyze bond-orientational order in a 2D image, restricted to a circular region
%
% Input:
%   I - 2D grayscale image (double or uint8)
%
% Outputs:
%   psi_global  - struct with fields psi2, psi4, psi6 (⟨|ψₙ|⟩ values)
%   order_ratio - scalar = ⟨|ψ₆|⟩ / ⟨|ψ₂|⟩
%   Npoints     - number of detected points used

    %% Step 1: Create mask
    [Ny, Nx] = size(I);
    % Parameters for elliptical mask
    a = Nx / 2.25;         % Semi-major axis (x-direction)
    b = Ny / 2.25;         % Semi-minor axis (y-direction)
    theta_deg = 0;     % Rotation angle in degrees
    theta = deg2rad(theta_deg);  % Convert to radians
    
    % Meshgrid coordinates
    [Ny, Nx] = size(I);
    [X, Y] = meshgrid(1:Nx, 1:Ny);
    cx = Nx / 2;
    cy = Ny / 2;
    
    % Shifted coordinates
    Xc = X - cx;
    Yc = Y - cy;
    
    % Rotate coordinates
    Xr =  cos(theta) * Xc + sin(theta) * Yc;
    Yr = -sin(theta) * Xc + cos(theta) * Yc;
    
    % Elliptical mask equation
    mask = (Xr / a).^2 + (Yr / b).^2 <= 1;
    
    % Apply mask to image
    I = double(I) .* mask;

    %% Step 2: Detect local maxima within mask
    I_smooth = imgaussfilt(I, 2);
    peaks = imregionalmax(I_smooth);
    [y, x] = find(peaks);
    Npoints = length(x);

    if Npoints < 6
        warning('Too few points (%d) to compute bond order meaningfully.', Npoints);
        psi_global = struct('psi2', NaN, 'psi4', NaN, 'psi6', NaN);
        order_ratio = NaN;
        return;
    end

    %% Step 3: Try Delaunay triangulation
    use_kNN = false;
    try
        tri = delaunay(x, y);
    catch
        warning('Delaunay triangulation failed - falling back to kNN neighbor search');
        use_kNN = true;
    end

    % 3) Find neighbors
    neighbors = cell(Npoints, 1);
    if ~use_kNN
        % Use Delaunay neighbors
        for i = 1:size(tri,1)
            for j = 1:3
                idx = tri(i,j);
                nbrs = tri(i,[1 2 3] ~= j);
                neighbors{idx} = unique([neighbors{idx}, nbrs]);
            end
        end
    else
        % Use kNN neighbors (e.g. k=6)
        k = min(6, Npoints-1);
        pts = [x, y];
        [idx, ~] = knnsearch(pts, pts, 'K', k+1); % +1 for self
        for j = 1:Npoints
            neighbors{j} = idx(j, 2:end);
        end
    end

    % 4) Compute bond order parameters
    psi2 = zeros(Npoints,1);
    psi4 = zeros(Npoints,1);
    psi6 = zeros(Npoints,1);

    for j = 1:Npoints
        nbrs = neighbors{j};
        nbrs(nbrs == j) = [];
        if isempty(nbrs)
            continue;
        end
        angles = atan2(y(nbrs) - y(j), x(nbrs) - x(j));
        psi2(j) = abs(mean(exp(1i * 2 * angles)));
        psi4(j) = abs(mean(exp(1i * 4 * angles)));
        psi6(j) = abs(mean(exp(1i * 6 * angles)));
    end

    psi_global.psi2 = mean(psi2);
    psi_global.psi4 = mean(psi4);
    psi_global.psi6 = mean(psi6);

    order_ratio = psi_global.psi6 / psi_global.psi2;

    % --- After neighbors are found and (x,y) extracted ---

    % 1) Angular distribution of all bonds pooled from all points
    all_angles = [];
    for j = 1:Npoints
        nbrs = neighbors{j};
        nbrs(nbrs == j) = [];
        if isempty(nbrs), continue; end
        angles = atan2(y(nbrs)-y(j), x(nbrs)-x(j));
        all_angles = [all_angles; angles(:)];
    end
    
    % Histogram of bond angles over [-pi, pi]
    nbins = 36; % 10 degree bins
    edges = linspace(-pi, pi, nbins+1);
    counts = histcounts(all_angles, edges, 'Normalization','probability');
    
    % Circular variance of bond angles
    R = abs(mean(exp(1i*all_angles)));
    circ_variance = 1 - R;  % 0 means angles clustered; 1 means uniform
    
    % 2) Positional correlation function g(r) estimate
    
    % Compute pairwise distances
    coords = [x y];
    D = pdist(coords);
    max_r = min([Nx, Ny])/2; % max radius for g(r)
    
    % Compute histogram of distances
    dr = 1; % bin width in pixels (adjust)
    r_edges = 0:dr:max_r;
    [counts_r, r_bins] = histcounts(D, r_edges);
    
    % Normalize g(r) by ideal gas distribution (circle annulus area)
    r_centers = r_edges(1:end-1) + dr/2;
    area_annulus = pi*((r_edges(2:end)).^2 - (r_edges(1:end-1)).^2);
    density = Npoints / (Nx*Ny);
    
    g_r = counts_r ./ (density * area_annulus * Npoints);

    figure(1);
    clf
    set(gcf,'Position',[500 100 1250 500])
    t = tiledlayout(1, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid

    nexttile
    imshow(I, []);
    hold on;
    scatter(x, y, 30, 'w', 'filled');
    hold off;
    colormap(Helper.Colormaps.plasma()); % Example colormap for original image plot
    cbar1 = colorbar;
    cbar1.Label.Interpreter = 'latex';
    xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14);
    ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14);
    title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14);
    axis square;

    % --- Plot angular distribution ---
    nexttile
    polarhistogram(all_angles, nbins, 'Normalization', 'probability');
    title(sprintf('Bond angle distribution (circ. variance = %.3f)', circ_variance));
    
    % --- Plot g(r) ---
    nexttile
    plot(r_centers, g_r, 'LineWidth', 2);
    xlabel('r (pixels)');
    ylabel('g(r)');
    title('Radial Pair Correlation Function g(r)');
    grid on;
    axis square


end