function [AveragePCD] = extractAveragePeakColumnDensity(folder_path, run_index, radius, minPeakHeight, SuppressPlotFlag)    
    set(0,'defaulttextInterpreter','latex')
    set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
    
    format long
    
    % 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);
        AveragePCD = NaN;
        return;
    end
    
    Params = Data.Params;
    Transf = Data.Transf; 
    Observ = Data.Observ; 
    
    if isgpuarray(Data.psi)
        psi = gather(Data.psi);
    else
        psi = Data.psi;
    end
    if isgpuarray(Data.Observ.residual)
        Observ.residual = gather(Data.Observ.residual); 
    else
        Observ.residual = Data.Observ.residual; 
    end
    
    % Axes scaling and coordinates 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 probability density |psi|^2
    n                       = abs(psi).^2;
    nxy                     = squeeze(trapz(n*dz,3));
    
    state                   = nxy';
    peaks                   = extractPeaks(state, radius, minPeakHeight);
    peakHeights             = state(peaks);
    AveragePCD              = mean(peakHeights);
    [row, col]              = find(peaks);

    if ~SuppressPlotFlag
        figure(1); 
        clf
        set(gcf,'Position', [100, 100, 900, 900])
        plotxy = pcolor(x,y,state);
        set(plotxy, 'EdgeColor', 'none');
        hold on;
        plot(x(col), y(row), 'w+', 'MarkerSize', 8, 'LineWidth', 1.5);
        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(['<CD> = ' num2str(AveragePCD)], ...
          'Interpreter', 'tex', 'FontSize', 16)
        set(gca, 'FontSize', 16);   % For tick labels only
    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