%% ========================================================================
%  MRI_bubble_pattern_analysis.m
%  ------------------------------------------------------------------------
%  Analysis of bubble patterns in vertically vibrated gas-fluidized beds
%  ========================================================================
%
%  This script accompanies the publication:
%
%      "Self-organization of bubbles into spatiotemporally regular
%       patterns in vertically vibrated 3D gas-fluidized granular beds"
%       Benjamin Rahimian1,2, Melis Özdemir1, Nick Hildebrandt3, Swantje Pietsch-Braune3, Stefan Benders1, Stefan Heinrich3, and Alexander Penn1*, 
%       Physical Review Letters, [year], DOI: [manuscript DOI]
%
%       1Institute of Process Imaging, Hamburg University of Technology, Hamburg, Germany
%       2Centre for Nature-Inspired Engineering and Department of Chemical Engineering, University College London, London, United Kingdom
%       3Institute of Solids Process Engineering and Particle Technology, Hamburg University of Technology, Hamburg, Germany
%
%  Associated data set:
%      TORE (TUHH Open Research), DOI: https://doi.org/10.15480/882.17457
%
%  DESCRIPTION
%    The script analyses MRI image series of a horizontal slice through a
%    vertically vibrated, gas-fluidized granular bed of poppy seeds. It contains two
%    independent analysis parts and the corresponding manuscript figures:
%
%      1) FFT Analysis
%         Computes the dominant bubbling frequency f_b of an image series
%         from the summed pixel intensity within an automatically
%         generated region-of-interest (ROI) mask. For each case, the
%         intensity time series and the normalized magnitude spectrum are
%         plotted. The section "Plot FFT results (Figure 2a)" reproduces
%         Fig. 2a of the manuscript from the tabulated results of all
%         measured cases.
%
%      2) Alternation Coefficient Analysis
%         Computes the alternation coefficient alpha of an image series by
%         comparing binarized bubble patterns one and two bubbling periods
%         apart using the Dice coefficient. The section "Plot
%         Alpha Coefficient Results" reproduces the corresponding figure
%         of the manuscript from the tabulated results of all measured
%         cases.
%
%  REQUIREMENTS
%    - MATLAB (developed and tested with R2023b)
%    - Image Processing Toolbox
%
%  DATA STRUCTURE
%    The data repository contains one example case (bed height H = 28 dp,
%    vibration frequency f_vib = 6 Hz) with the following layout:
% 
%        <dataRoot>/                        e.g. .../MRI_Scans
%            H28dp_f6Hz/
%                Horizontal_Slice/
%                    H28dp_f6Hz_horizontal_frames/
%                        <numbered PNG frames of the image series>
%
%    The PNG frames are grayscale images numbered consecutively in their
%    filename; they are sorted numerically by the script.
%
%  USAGE
%    1) Download the data set from the repository (DOI above).
%    2) In the "Read Data" section of each analysis part, set 'dataRoot'
%       to the location of the MRI_Scans folder on your machine. If your
%       PNG frames sit in a further subfolder, extend 'framePath'
%       accordingly.
%    3) Run the script section by section ("Run Section" / Ctrl+Enter).
%       Note: both analysis parts begin with 'clear all' and are fully
%       independent of each other.
%    4) The processing loops accept an arbitrary number of cases. To
%       analyse additional image series, add further entries to the
%       'cases' struct array in the "Read Data" sections (see comments
%       there). The two figure sections that reproduce the manuscript
%       figures use hard-coded tabulated results of ALL measured cases,
%       since the repository contains only one example image series.
%
%  CITATION
%    If you use this code and/or the associated data, please cite the
%    publication given above.
%
%  AUTHOR
%    Benjamin Rahimian (benjamin.rahimian@gmail.com)
%
%  VERSION
%    v1.0 - 09.07.2026
%% ========================================================================

%% FFT Analysis
clear all
clc
close all

%% Read Data (FFT analysis)

% ---- User settings ------------------------------------------------------
% Root folder of the downloaded data set. Adjust to your local path.
dataRoot = 'MRI_Scans';

% Cases to process. Each case needs:
%   name      - identifier used in printed output and result fields
%   framePath - folder containing the PNG image series of that case
%   f_vib     - vibration frequency of the bed in Hz
% The loop below works for any number of cases (including a single one);
% add further entries (cases(2), cases(3), ...) to analyse more series.
cases = struct('name', {}, 'framePath', {}, 'f_vib', {});

cases(1).name      = 'H28dp_f6Hz';
cases(1).framePath = fullfile(dataRoot, 'H28dp_f6Hz', 'Horizontal_Slice', ...
                              'H28dp_f6Hz_horizontal_frames');
cases(1).f_vib     = 6;

dt = 0.03; % Temporal resolution of the image series in seconds
% -------------------------------------------------------------------------

% Process each case
for caseIdx = 1:numel(cases)
    caseName    = cases(caseIdx).name;
    frameFolder = cases(caseIdx).framePath;
    f_vib       = cases(caseIdx).f_vib;

    fprintf('\n=== Processing %s (f_vib = %g Hz) ===\n', caseName, f_vib);

    if ~isfolder(frameFolder)
        fprintf('Folder not found: %s. Skipping...\n', frameFolder);
        continue;
    end

    % Get PNG files of the image series
    fileInfo = dir(fullfile(frameFolder, '*.png'));
    if isempty(fileInfo)
        fprintf('No PNGs found in %s. Skipping...\n', frameFolder);
        continue;
    end

    % Sort frames numerically by the number contained in the filename
    fileNumbers = arrayfun(@(f) str2double(regexp(f.name, '\d+', 'match', 'once')), fileInfo);
    [~, sortIdx] = sort(fileNumbers);
    fileInfo = fileInfo(sortIdx);

    numFrames = length(fileInfo);

    % === Build ROI mask from SUM of all frame intensities ===
    % Sum in double to avoid uint8 saturation, then rescale to [0, 255]
    sumImg = [];
    for k = 1:numFrames
        img = imread(fullfile(frameFolder, fileInfo(k).name));
        if size(img, 3) == 3
            img = rgb2gray(img);
        end
        if isempty(sumImg)
            sumImg = double(img);
        else
            sumImg = sumImg + double(img);
        end
    end
    sumImgNorm = sumImg - min(sumImg(:));
    if max(sumImgNorm(:)) > 0
        sumImgNorm = sumImgNorm / max(sumImgNorm(:));
    end
    roiImage = uint8(round(255 * sumImgNorm));
    roiImage(roiImage < 50) = 0;
    roiFiltered = imgaussfilt(roiImage, 0.5);
    roiMask = imbinarize(roiFiltered, "adaptive", ...
        "ForegroundPolarity", "bright", ...
        "Sensitivity", 1);
    roiMask = bwareaopen(roiMask, 50);
    roiMask = imfill(roiMask, 'holes');

    % === Process all frames - summed pixel intensity within the ROI ===
    summedIntensitySignal = zeros(1, numFrames);

    fprintf('Processing %d frames...\n', numFrames);

    for i = 1:numFrames
        img = imread(fullfile(frameFolder, fileInfo(i).name));

        if size(img, 3) == 3
            img = rgb2gray(img);
        end

        % Raw summed intensity within mask
        maskedPixels = double(img(roiMask));
        summedIntensitySignal(i) = sum(maskedPixels);

    end

    % === FFT Analysis ===
    % Remove DC component
    signal = summedIntensitySignal - mean(summedIntensitySignal);

    % FFT
    Y = fft(signal);
    N = length(signal);

    % Frequency vector (positive frequencies only)
    fs = 1/dt;
    f = (0:N-1) * fs / N;
    f = f(1:floor(N/2)+1);

    % Magnitude spectrum
    magnitude = abs(Y(1:floor(N/2)+1));

    % Find dominant frequency (skip DC)
    [~, maxIdx] = max(magnitude(2:end));
    dominantFreq = f(maxIdx + 1);

    % Results
    fprintf('\n%s: Dominant frequency = %.4f Hz\n', caseName, dominantFreq);
    fprintf('Period = %.2f seconds (%.1f frames)\n', 1/dominantFreq, 1/(dominantFreq*dt));

    % === Plotting ===
    % Normalize x-axis by vibration frequency
    f_norm = f / f_vib;

    % Normalize magnitude so that peak is 1
    magnitude_normalized_plot = magnitude / max(magnitude);

    % Style parameters
    lineWidth = 1.2;
    fontSize = 16;

    % --- Plot time series ---
    figure('Color', 'w', 'Units', 'normalized', 'Position', [0.1, 0.5, 0.4, 0.25]);
    % Time axis normalized by the total scan duration t_scan = numFrames*dt
    timeVec = ((0:N-1)*dt) / (numFrames*dt);
    plot(timeVec, summedIntensitySignal, 'k', 'LineWidth', lineWidth);
    xlabel('t/t_{scan}', 'FontSize', fontSize, 'Color', 'k');
    ylabel('I(n)', 'FontSize', fontSize, 'Color', 'k');
    set(gca, ...
        'Box', 'on', ...
        'FontSize', fontSize, ...
        'LineWidth', 1, ...
        'XColor', 'k', ...
        'YColor', 'k');
    axis tight;

    % --- Plot normalized FFT ---
    figure('Color', 'w', 'Units', 'normalized', 'Position', [0.5, 0.5, 0.4, 0.25]);
    plot(f_norm(2:end), magnitude_normalized_plot(2:end), 'k', 'LineWidth', lineWidth);
    xlabel('f/f_0', 'FontSize', fontSize, 'Color', 'k');
    ylabel('|X(m)|(max|X(m)|)^{-1}', 'FontSize', fontSize, 'Color', 'k');
    set(gca, ...
        'Box', 'on', ...
        'FontSize', fontSize, ...
        'LineWidth', 1, ...
        'XColor', 'k', ...
        'YColor', 'k');
    % Optional: reduce number of ticks
    xticks(0:1:floor(max(f_norm)));
    yticks(0:0.5:1);

    axis tight;
    hold on;

    % Reference lines at the subharmonic, fundamental and harmonics
    xline(0.5, '--', 'Color', [0.4 0.4 0.4], 'LineWidth', lineWidth);  % dashed at 0.5 f0
    xline(1,   '-',  'Color', [0.4 0.4 0.4], 'LineWidth', lineWidth);  % solid at f0
    xline(2,   '--', 'Color', [0.4 0.4 0.4], 'LineWidth', lineWidth);  % dashed at 2 f0
    xline(3,   '--', 'Color', [0.4 0.4 0.4], 'LineWidth', lineWidth);  % dashed at 3 f0

    % Mark dominant frequency peak
    dom_freq_normalized = dominantFreq / f_vib;
    dom_peak_idx = find(abs(f_norm - dom_freq_normalized) == min(abs(f_norm - dom_freq_normalized)), 1);
    if ~isempty(dom_peak_idx)
        plot(dom_freq_normalized, magnitude_normalized_plot(dom_peak_idx), 'r*', ...
             'MarkerSize', 8,'LineWidth', 2, 'MarkerFaceColor', 'none');
    end
    ylim([0 1.1]);

    % Store results for this case
    % The summary figure below uses hard-coded tabulated values of all
    % measured cases instead. Kept for verification before removal.
    results.(caseName).dominantFreq = dominantFreq;
    results.(caseName).f_vib = f_vib;
    results.(caseName).ratio = dominantFreq / f_vib;
    results.(caseName).numFrames = numFrames;
    results.(caseName).totalTime = numFrames * dt;

    fprintf('Frequency ratio (f_dom/f_vib): %.3f\n', dominantFreq / f_vib);
end

%% Plot FFT results (Figure 2a)
% Reproduces Fig. 2a of the manuscript from the tabulated dominant
% bubbling frequencies of ALL measured cases (bed heights H = 28 dp and
% H = 47 dp). These values are hard-coded because the data repository
% contains only one example image series.

% --- Data ---
% H=28dp data
f_vib_H28 = [0, 3, 4, 5, 6, 7, 8, 9];
f_dom_H28 = [9.5556, 5.8333, 3.8889, 4.8333, 5.8333, 3.3889, 3.8889, 4.4444];

% H=47dp data
f_vib_H47 = [0, 3, 4, 5, 6, 7, 8, 9];
f_dom_H47 = [4.6667, 2.8889, 3.8889, 4.8333, 5.8333, 3.3889, 3.5000, 3.5556];

% --- Color palette ---
color_H28 = [0.25, 0.55, 0.96];  % Darker blue
color_H47 = [0.96, 0.26, 0.31]; % Darker red
grey_line = [0.45, 0.45, 0.45];          % Grey for reference lines

% === Main Figure ===
mainFig = figure('Color', 'w');
hold on;

% --- Plot reference lines first (behind data) ---
f_range = [0, 10];
% y = x line
plot(f_range, f_range, '--', 'Color', grey_line, 'LineWidth', 2, ...
    'HandleVisibility', 'off');
% y = 2x line
plot(f_range, 2*f_range, '--', 'Color', grey_line, 'LineWidth', 2, ...
    'HandleVisibility', 'off');
% y = 0.5x line
plot(f_range, 0.5*f_range, '--', 'Color', grey_line, 'LineWidth', 2, ...
    'HandleVisibility', 'off');

% --- Plot data points (excluding the unvibrated reference f_vib = 0) ---
% H=28dp
h1 = plot(f_vib_H28(2:end), f_dom_H28(2:end), ...
    'Marker', '+', ...
    'LineStyle', 'none', ...
    'MarkerSize', 15, ...
    'MarkerEdgeColor', color_H28, ...
    'MarkerFaceColor', 'none', ...
    'LineWidth', 2, ...
    'DisplayName', 'H = 28 dp');

% H=47dp
h2 = plot(f_vib_H47(2:end), f_dom_H47(2:end), ...
    'Marker', 'x', ...
    'LineStyle', 'none', ...
    'MarkerSize', 15, ...
    'MarkerEdgeColor', color_H47, ...
    'MarkerFaceColor', 'none', ...
    'LineWidth', 2, ...
    'DisplayName', 'H = 47 dp');

% --- Axes labels and formatting ---
xlabel('f_{vib} / Hz', 'FontSize', 20, 'FontName', 'Times New Roman');
ylabel('f_b / Hz', 'FontSize', 20, 'FontName', 'Times New Roman');
xlim([2.5, 9.5]);
ylim([2.5, 6.5]);
set(gca, 'LineWidth', 1);
set(gca, 'FontSize', 20, 'FontName', 'Times New Roman');
box on;
grid off;

% --- Set equal axis scaling ---
axis equal;
xlim([2.5, 9.5]);
ylim([2.5, 6.5]);

% --- Legend ---
lgd = legend([h1, h2], 'Location', 'northeast');
set(lgd, 'FontSize', 18, 'FontName', 'Times New Roman', 'Box', 'off');

hold off;

%% Alternation coefficient Analysis
clear all
clc
close all

%% Read Data (alternation coefficient analysis)

% ---- User settings ------------------------------------------------------
% Root folder of the downloaded data set. Adjust to your local path.
dataRoot = 'MRI_Scans';

% Cases to process. Each case needs:
%   name      - identifier used in printed output
%   framePath - folder containing the PNG image series of that case
%   f_dom     - dominant bubbling frequency of that case in Hz, as
%               obtained from the FFT analysis above
% The loop below works for any number of cases (including a single one);
% add further entries (cases(2), cases(3), ...) to analyse more series.
cases = struct('name', {}, 'framePath', {}, 'f_dom', {});

cases(1).name      = 'H28dp_f6Hz';
cases(1).framePath = fullfile(dataRoot, 'H28dp_f6Hz', 'Horizontal_Slice', ...
                              'H28dp_f6Hz_horizontal_frames');
cases(1).f_dom     = 5.8333;

dt = 0.03; % Temporal resolution of the image series in seconds
% -------------------------------------------------------------------------

% Result storage: [case index, mean alpha, std of alpha]
results = zeros(numel(cases), 3);

for idx = 1:numel(cases)
    caseName    = cases(idx).name;
    frameFolder = cases(idx).framePath;

    if ~isfolder(frameFolder)
        warning('Folder not found: %s', frameFolder);
        continue;
    end

    % Get PNG files of the image series
    fileInfo = dir(fullfile(frameFolder, '*.png'));
    if isempty(fileInfo)
        warning('No PNGs found in %s', frameFolder);
        continue;
    end

    % Sort frames numerically by the number contained in the filename
    fileNumbers = arrayfun(@(f) str2double(regexp(f.name, '\d+', 'match', 'once')), fileInfo);
    [~, sortIdx] = sort(fileNumbers);
    fileInfo = fileInfo(sortIdx);
    numFrames = length(fileInfo);

    % === Build ROI mask from SUM of all frame intensities ===
    % Sum in double to avoid uint8 saturation, then rescale to [0, 255]
    sumImg = [];
    for k = 1:numFrames
        img = imread(fullfile(frameFolder, fileInfo(k).name));
        if size(img, 3) == 3
            img = rgb2gray(img);
        end
        if isempty(sumImg)
            sumImg = double(img);
        else
            sumImg = sumImg + double(img);
        end
    end

    sumImgNorm = sumImg - min(sumImg(:));
    if max(sumImgNorm(:)) > 0
        sumImgNorm = sumImgNorm / max(sumImgNorm(:));
    end
    roiImage = uint8(round(255 * sumImgNorm));

    roiImage(roiImage < 50) = 0;
    roiFiltered = imgaussfilt(roiImage, 0.5);
    roiMask = imbinarize(roiFiltered, "adaptive", ...
        "ForegroundPolarity", "bright", ...
        "Sensitivity", 1);

    roiMask = bwareaopen(roiMask, 50);
    roiMask = imfill(roiMask, 'holes');

    % shrink mask by 1-pixel layer
    roiMask = imerode(roiMask, strel('disk', 1));

    % === Load all frames ===
    % (raw frames are stored here; binarization happens in the next step)
    extractedROI = cell(numFrames, 1);
    for i = 1:numFrames
        img = imread(fullfile(frameFolder, fileInfo(i).name));
        extractedROI{i} = img;
    end

    % === Binary Processing ===
    % Threshold and binarize each frame, then keep only the dark (bubble)
    % regions that lie fully inside the ROI (i.e. not connected to the
    % background outside the bed).
    BWfinal = cell(numFrames, 1);

    for i = 1:numFrames
        J = extractedROI{i};
        J(J < 70) = 0;
        J2 = imgaussfilt(J, 0.5);
        BW = imbinarize(J2, "adaptive", ...
                    "ForegroundPolarity", "bright", ...
                    "Sensitivity", 0.8);

        inverted = ~BW;
        blackInMask = inverted & roiMask;
        backgroundMask = ~roiMask;
        backgroundConnected = imreconstruct(backgroundMask, blackInMask);
        cleanBlackSpots = blackInMask & ~backgroundConnected;

        BWfinal{i} = bwareaopen(cleanBlackSpots, 8);
    end

    % === DICE Score Analysis with 3-way Comparison ===
    % Compare binarized bubble patterns one and two dominant bubbling
    % periods apart; alpha is the mean difference of the two Dice scores.
    f_dom = cases(idx).f_dom;                 % dominant bubbling frequency in Hz
    frameRate = 1/dt;
    exactFramesPerPeriod = frameRate / f_dom; % frames per bubbling period

    % Define the two comparison distances
    framesPer1Period = round(1*exactFramesPerPeriod);  % 1 period
    framesPer2Period = round(2*exactFramesPerPeriod);  % 2 periods

    % Split the series into three sections to estimate a standard deviation
    sectionLength = floor(numFrames / 3);
    sections = {1:sectionLength, sectionLength+1:2*sectionLength, 2*sectionLength+1:numFrames};

    sectionMeans = zeros(1, 3);
    for s = 1:3
        indices = sections{s};
        diceMetrics = [];

        % Calculate for frames that allow both comparisons
        maxFrameIndex = min(indices(end), numFrames - framesPer2Period);

        for i = indices(1):maxFrameIndex
            A = logical(BWfinal{i});
            B1 = logical(BWfinal{i + framesPer1Period}); % 1 period apart
            B2 = logical(BWfinal{i + framesPer2Period}); % 2 periods apart

            totalWhiteA = sum(A(:));
            totalWhiteB1 = sum(B1(:));
            totalWhiteB2 = sum(B2(:));

            % Check if all comparisons are valid (at least one frame has objects)
            valid1Period = (totalWhiteA > 0 || totalWhiteB1 > 0);
            valid2Period = (totalWhiteA > 0 || totalWhiteB2 > 0);

            if valid1Period && valid2Period
                % Calculate DICE for 1 period apart
                if totalWhiteA == 0 || totalWhiteB1 == 0
                    dice1Period = 0;
                else
                    intersection1 = sum(A(:) & B1(:));
                    total1 = totalWhiteA + totalWhiteB1;
                    dice1Period = 2 * intersection1 / total1;
                end

                % Calculate DICE for 2 periods apart
                if totalWhiteA == 0 || totalWhiteB2 == 0
                    dice2Period = 0;
                else
                    intersection2 = sum(A(:) & B2(:));
                    total2 = totalWhiteA + totalWhiteB2;
                    dice2Period = 2 * intersection2 / total2;
                end

                diceMetric = (dice2Period - dice1Period);
                diceMetrics(end+1) = diceMetric;
            end
        end

        if ~isempty(diceMetrics)
            sectionMeans(s) = mean(diceMetrics);
        else
            sectionMeans(s) = NaN;
        end
    end

    % Final stats
    validMeans = sectionMeans(~isnan(sectionMeans));
    meanDiceMetric = mean(validMeans);
    stdDiceMetric = std(validMeans);

    % Store result
    results(idx, :) = [idx, meanDiceMetric, stdDiceMetric];

    % Print
    fprintf('%s: Mean alpha = %.4f | Std = %.4f | ROI from sum of %d frames\n', ...
        caseName, meanDiceMetric, stdDiceMetric, numFrames);
end

%% Plot Alpha Coefficent Results
% Reproduces the alpha-coefficient figure of the manuscript from the
% tabulated results of ALL measured cases (bed heights H = 28 dp and
% H = 47 dp). These values are hard-coded because the data repository
% contains only one example image series.

% Data matrix: [Mean_Metric, Std, Frequency, Group, Pattern]

data = [
    % Group 1: H=33mm=28dp (Blue)
    0.0588, 0.0101, 3, 1, 0;  % f=3, no pattern
    0.5799, 0.0522, 4, 1, 2;  % f=4, 2-2 pattern
    0.5497, 0.0233, 5, 1, 1;  % f=5, 1-1:f/2 pattern
    0.5400, 0.0068, 6, 1, 2;  % f=6, 2-2 pattern
    0.1669, 0.0314, 7, 1, 0;  % f=7, no pattern
    0.4609, 0.0346, 8, 1, 3;  % f=8, 1-1:f/4 pattern
    0.1700, 0.0572, 9, 1, 0;  % f=9, no pattern

    % Group 2: H=55mm=47dp (Red)
    0.0646, 0.0091, 3, 2, 0;  % f=3, no pattern
    0.6850, 0.0173, 4, 2, 1;  % f=4, 1-1:f/2 pattern
    0.1391, 0.0809, 5, 2, 0;  % f=5, no pattern
    0.0175, 0.0145, 6, 2, 0;  % f=6, no pattern
    0.5019, 0.0249, 7, 2, 3;  % f=7, 1-1:f/4 pattern
    0.0097, 0.0277, 8, 2, 0;  % f=8, no pattern
    0.0258, 0.0218, 9, 2, 0;  % f=9, no pattern
];

% --- Extract data ---
y = data(:,1);           % Mean Metric
std_dev = data(:,2);     % Standard deviation
x = data(:,3);           % Frequency
group = data(:,4);       % Group number (1=H=28dp, 2=H=47dp)
pattern = data(:,5);     % Pattern type (0=none, 1=1-1:f/2, 2=2-2, 3=1-1:f/4)

% --- Color palette ---
color1 = [0.25, 0.55, 0.96];   % blue
color2 = [0.96, 0.26, 0.31];   % red

% --- Assign fill colors based on group ---
fillColors = cell(length(y), 1);
for i = 1:length(y)
    if group(i) == 1
        fillColors{i} = color1;
    else
        fillColors{i} = color2;
    end
end

% --- Marker shapes by pattern ---
marker_map = containers.Map([0, 1, 2, 3], {'o', 's', 'd', '^'});

% --- Assign markers based on pattern ---
markers = cell(length(y), 1);
for i = 1:length(y)
    markers{i} = marker_map(pattern(i));
end

% --- Display names for legend ---
customNames = {
    'H = 28 dp', ...
    'H = 47 dp'
};

% Pattern names for legend
patternNames = {
    'No pattern', ...
    '1-1:f/2 Pattern', ...
    '2-2 Pattern', ...
    '1-1:f/4 Pattern'
};

% === Main Figure ===
mainFig = figure('Color', 'w');
hold on;

% y-axis limits
y_min = -0.05;
y_max = 0.75;

% Plot markers only
for i = 1:length(y)
    plot(x(i), y(i), ...
        'Marker', markers{i}, ...
        'LineStyle', 'none', ...
        'MarkerSize', 15, ...
        'MarkerEdgeColor', 'k', ...
        'MarkerFaceColor', fillColors{i}, ...
        'LineWidth', 1, ...
        'HandleVisibility', 'off');
end

% Plot error bars on top
for i = 1:length(y)
    errorbar(x(i), y(i), std_dev(i), ...
        'LineStyle', 'none', ...
        'Color', 'k', ...
        'LineWidth', 1, ...
        'HandleVisibility', 'off', ...
        'Marker', 'none');
end

% Axes labels and formatting
xlabel('f_{vib} / Hz', 'FontSize', 18, 'FontName', 'Times New Roman');
ylabel('\alpha / -', 'FontSize', 18, 'FontName', 'Times New Roman');
set(gca, 'XTick', unique(x));
ylim([y_min y_max]);
xlim([2.5 9.5]);
set(gca, 'LineWidth', 1);
set(gca, 'FontSize', 18, 'FontName', 'Times New Roman');
box on;
grid off;
set(gca, 'Layer', 'top');

% === Create SEPARATE legend figure ===

% Create a new figure for the legend only
legendFig = figure('Color', 'w');
hold on;

% Create dummy plots with the same properties
h1 = plot(nan, nan, ...
    'Marker', 'o', ...
    'LineStyle', 'none', ...
    'MarkerSize', 15, ...
    'Color', 'k', ...
    'LineWidth', 1, ...
    'MarkerEdgeColor', color1, ...
    'MarkerFaceColor', color1, ...
    'DisplayName', customNames{1});

h2 = plot(nan, nan, ...
    'Marker', 'o', ...
    'LineStyle', 'none', ...
    'MarkerSize', 15, ...
    'Color', 'k', ...
    'LineWidth', 1, ...
    'MarkerEdgeColor', color2, ...
    'MarkerFaceColor', color2, ...
    'DisplayName', customNames{2});

h3 = plot(nan, nan, ...
    'Marker', 'o', ...
    'LineStyle', 'none', ...
    'MarkerSize', 15, ...
    'Color', 'k', ...
    'LineWidth', 1, ...
    'MarkerEdgeColor', 'k', ...
    'MarkerFaceColor', "none", ...
    'DisplayName', patternNames{1});

h4 = plot(nan, nan, ...
    'Marker', 's', ...
    'LineStyle', 'none', ...
    'MarkerSize', 15, ...
    'Color', 'k', ...
    'LineWidth', 1, ...
    'MarkerEdgeColor', 'k', ...
    'MarkerFaceColor', "none", ...
    'DisplayName', patternNames{2});

h5 = plot(nan, nan, ...
    'Marker', 'd', ...
    'LineStyle', 'none', ...
    'MarkerSize', 15, ...
    'Color', 'k', ...
    'LineWidth', 1, ...
    'MarkerEdgeColor', 'k', ...
    'MarkerFaceColor', "none", ...
    'DisplayName', patternNames{3});

h6 = plot(nan, nan, ...
    'Marker', '^', ...
    'LineStyle', 'none', ...
    'MarkerSize', 15, ...
    'Color', 'k', ...
    'LineWidth', 1, ...
    'MarkerEdgeColor', 'k', ...
    'MarkerFaceColor', "none", ...
    'DisplayName', patternNames{4});

% Create the legend (use 'north' or 'best' instead of 'center')
lgd = legend([h1, h2, h3, h4, h5, h6], 'Location', 'north', 'NumColumns', 1);
set(lgd, 'FontSize', 18, 'FontName', 'Times New Roman', 'Box', 'off');
axis off;
set(legendFig, 'Position', [100, 100, 400, 350]);
hold off;