function [ordered_od_imgs, ordered_scan_parameter_values, ordered_file_list] = collectODImages(options)
%% Applies cropping, background subtraction, and optional fringe removal, optional unshuffling on OD image dataset
% Automatically reuses in-memory full dataset if available; 
% otherwise, reads and processes raw HDF5 data.
%
% Inputs:
%   options - structure containing processing options:
%       .folderPath           : path to raw HDF5 files
%       .saveDirectory        : path to save cache (if needed)
%       .cam, .angle          : camera selection and rotation angle
%       .ImagingMode, .PulseDuration : imaging parameters
%       .scan_parameter       : name of scan parameter
%       .center, .span        : cropping settings
%       .fraction             : background subtraction fraction
%       .removeFringes        : logical flag for fringe removal
%       .skipUnshuffling      : logical flag to skip unshuffling
%       .scan_reference_values: reference values for unshuffling
%
% Outputs:
%   ordered_od_imgs              : cell array of processed OD images (ordered)
%   ordered_scan_parameter_values: vector of scan parameter values (ordered)
%   ordered_file_list            : cell array of file names (ordered)
    
    % --- Early exit if processed data already exist ---
    if evalin('base', 'exist(''od_imgs'',''var'') && exist(''scan_parameter_values'',''var'') && exist(''file_list'',''var'')')
        fprintf('\nReusing processed OD images, scan parameters, and file list from memory.\n');
        ordered_od_imgs               = evalin('base','od_imgs');
        ordered_scan_parameter_values = evalin('base','scan_parameter_values');
        ordered_file_list             = evalin('base','file_list');
        return; % ✅ skip rest of the function
    end

    % --- Check if the full OD dataset and scan parameters exist in workspace ---
    fullDataExists = evalin('base', 'exist(''full_od_imgs'', ''var'')') && ...
                     evalin('base', 'exist(''full_bkg_imgs'', ''var'')') && ...
                     evalin('base', 'exist(''raw_scan_parameter_values'', ''var'')') && ...
                     evalin('base', 'exist(''raw_file_list'', ''var'')');

    if fullDataExists
        % Both required datasets exist, use them directly
        fprintf('\nReusing full OD image dataset and scan parameters from memory.\n');
        full_od_imgs              = evalin('base', 'full_od_imgs');
        full_bkg_imgs             = evalin('base', 'full_bkg_imgs');
        raw_scan_parameter_values = evalin('base', 'raw_scan_parameter_values');
        raw_file_list             = evalin('base', 'raw_file_list');
    else
        % Either dataset is missing, process raw HDF5 files completely
        fprintf('\nFull OD image dataset or scan parameters not found in memory.\n');
        [full_od_imgs, full_bkg_imgs, raw_scan_parameter_values, raw_file_list] = Helper.processRawData(options);

        % Save raw full dataset for reuse
        assignin('base', 'full_od_imgs', full_od_imgs);
        assignin('base', 'full_bkg_imgs', full_bkg_imgs);
        assignin('base', 'raw_scan_parameter_values', raw_scan_parameter_values);
        assignin('base', 'raw_file_list', raw_file_list);
        fprintf('\nCompleted computing OD images. Stored in workspace for reuse.\n');
    end

    nFiles = size(full_od_imgs, 3);

    % --- Preallocate arrays for processed images ---
    absimages = zeros(options.span(1)+1, options.span(2)+1, nFiles, 'single');
    refimages = zeros(options.span(1)+1, options.span(2)+1, nFiles, 'single');

    % --- Process each image: crop and subtract background ---
    for k = 1:nFiles
        od_img  = full_od_imgs(:,:,k);   % original full OD image, never modified
        bkg_img = full_bkg_imgs(:,:,k);  % original full background image, never modified
        if any(isnan(od_img(:)))
            absimages(:,:,k) = nan(options.span(1)+1, options.span(2)+1, 'single');
            continue
        end
        if any(isnan(bkg_img(:)))
            refimages(:,:,k) = nan(options.span(1)+1, options.span(2)+1, 'single');
            continue
        end

        % Crop image around the region of interest
        cropped_absimage   = Helper.cropODImage(od_img, options.center, options.span);
        cropped_refimage   = Helper.cropODImage(bkg_img, options.center, options.span);

        % Subtract background offset based on fraction
        processed_absimage = Helper.subtractBackgroundOffset(cropped_absimage, options.fraction);
        processed_refimage = Helper.subtractBackgroundOffset(cropped_refimage, options.fraction);

        % Store processed image (transpose to match orientation)
        absimages(:,:,k) = processed_absimage';
        refimages(:,:,k) = processed_refimage';
    end

    % --- Optional fringe removal ---
    if isfield(options, 'removeFringes') && options.removeFringes
        fprintf('\nApplying fringe removal to processed images...\n');
        optrefimages             = Helper.removeFringesInImage(absimages, refimages);
        absimages_fringe_removed = absimages - optrefimages;
        processed_od_imgs        = arrayfun(@(i) absimages_fringe_removed(:,:,i), 1:nFiles, 'UniformOutput', false);
        fprintf('\nFringe removal completed.\n');
    else
        processed_od_imgs        = arrayfun(@(i) absimages(:,:,i), 1:nFiles, 'UniformOutput', false);
    end

    % --- Optional unshuffling based on scan reference values ---
    if isfield(options, 'skipUnshuffling') && ~options.skipUnshuffling
        fprintf('\nReordering images according to scan parameter reference values...\n');

        n_values = length(options.scan_reference_values);
        n_total  = length(raw_scan_parameter_values);
        n_reps   = n_total / n_values;

        ordered_scan_parameter_values = zeros(1, n_total);
        ordered_od_imgs     = cell(1, n_total);
        ordered_file_list   = cell(1, n_total);
        counter = 1;

        temp_scan_values = raw_scan_parameter_values;
        temp_od_imgs     = processed_od_imgs;
        temp_file_list   = raw_file_list;

        for rep = 1:n_reps
            for val = options.scan_reference_values
                idx = find(temp_scan_values == val, 1, 'first');
                if isempty(idx), continue; end
                ordered_scan_parameter_values(counter) = temp_scan_values(idx);
                ordered_od_imgs{counter} = temp_od_imgs{idx};
                ordered_file_list{counter} = temp_file_list{idx};
                temp_scan_values(idx) = NaN;
                temp_od_imgs{idx} = [];
                temp_file_list{idx} = [];
                counter = counter + 1;
            end
        end
        fprintf('\nImage reordering completed.\n');
    else
        % No unshuffling: keep original order
        ordered_od_imgs = processed_od_imgs;
        ordered_scan_parameter_values = raw_scan_parameter_values;
        ordered_file_list             = raw_file_list;
    end
    
    % --- Save processed dataset for reuse ---
    assignin('base', 'od_imgs', ordered_od_imgs);
    assignin('base', 'scan_parameter_values', ordered_scan_parameter_values);
    assignin('base', 'file_list', ordered_file_list);

    fprintf('\nOD image dataset ready for further analysis.\n');
end