function [od_imgs, scan_parameter_values, 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:
%   od_imgs              : cell array of processed OD images
%   scan_parameter_values: vector of scan parameter values
%   file_list            : cell array of file names
    
    % --- Early exit if processed data already exist AND options match ---
    reuseVarsExist = evalin('base', ...
        'exist(''od_imgs'',''var'') && exist(''scan_parameter_values'',''var'') && exist(''file_list'',''var'') && exist(''prior_options'',''var'')');

    if reuseVarsExist
        prior_options = evalin('base','prior_options');

        % Define which fields are critical for reuse
        critical_fields = {'folderPath','cam','angle','ImagingMode','PulseDuration','center','span','fraction','removeFringes','skipUnshuffling','scan_reference_values'};

        if ~haveOptionsChanged(options, prior_options, critical_fields)
            fprintf('\nReusing processed OD images, scan parameters, and file list from memory (options unchanged).\n');
            od_imgs               = evalin('base','od_imgs');
            scan_parameter_values = evalin('base','scan_parameter_values');
            file_list             = evalin('base','file_list');
            return; % ✅ skip rest of the function
        else
            fprintf('\nProcessed-data-related options changed. Reprocessing full OD image dataset...\n');
        end
    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, check if raw-data options changed
        prior_options = evalin('base','prior_options');

        % Define critical fields that affect raw-data computation
        critical_raw_fields = {'folderPath','cam','angle','ImagingMode','PulseDuration'};

        if ~haveOptionsChanged(options, prior_options, critical_raw_fields)
            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
            fprintf('\nRaw-data-related options changed. Recomputing full OD image dataset...\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 recomputing OD images. Stored in workspace for reuse.\n');
        end
    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 and will be 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
        full_od_img  = full_od_imgs(:,:,k);   % original full OD image, never modified
        full_bkg_img = full_bkg_imgs(:,:,k);  % original full background image, never modified
        if any(isnan(full_od_img(:)))
            absimages(:,:,k) = nan(options.span(1)+1, options.span(2)+1, 'single');
            continue
        end
        if any(isnan(full_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(full_od_img, options.center, options.span);
        cropped_refimage   = Helper.cropODImage(full_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;
        od_imgs                  = arrayfun(@(i) absimages_fringe_removed(:,:,i), 1:nFiles, 'UniformOutput', false);
        fprintf('\nFringe removal completed.\n');
    else
        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     = 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
        od_imgs                       = ordered_od_imgs;
        scan_parameter_values         = ordered_scan_parameter_values;
        file_list                     = ordered_file_list;
        fprintf('\nImage reordering completed.\n');
    else
        % No unshuffling: keep original order
        scan_parameter_values         = raw_scan_parameter_values;
        file_list                     = raw_file_list;
    end
    
    % --- Save processed dataset and options for reuse ---
    assignin('base', 'od_imgs', od_imgs);
    assignin('base', 'scan_parameter_values', scan_parameter_values);
    assignin('base', 'file_list', file_list);
    assignin('base', 'prior_options', options);
    
    % --- Optionally save OD images as figures ---
    if ~options.skipSaveFigures
        odFolder = fullfile(saveDirectory, "Results", "ODImages");
        if ~exist(odFolder, 'dir')
            mkdir(odFolder);
        end
    
        for k = 1:length(od_imgs)
            img = od_imgs{k};
            fileName = fullfile(odFolder, sprintf('OD_img_%03d.png', k));
            imwrite(mat2gray(img), fileName);
        end
    end

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

end

% --- Local helper function to compare options ---
function changed = haveOptionsChanged(options, prior_options, critical_fields)
    changed = false;
    for f = critical_fields
        fname = f{1};
        if isfield(options, fname) && isfield(prior_options, fname)
            if ~isequal(options.(fname), prior_options.(fname))
                changed = true; return
            end
        elseif xor(isfield(options, fname), isfield(prior_options, fname))
            changed = true; return
        end
    end
end