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