function [full_od_imgs, full_bkg_imgs, raw_scan_parameter_values, raw_file_list] = processRawData(options)
%% Reads HDF5 files, computes OD images, supports disk-backed storage in blocks

    fprintf('\n[INFO] Processing raw data files at %s ...\n', options.folderPath);

    groupList = ["/images/MOT_3D_Camera/in_situ_absorption", ...
                 "/images/ODT_1_Axis_Camera/in_situ_absorption", ...
                 "/images/ODT_2_Axis_Camera/in_situ_absorption", ...
                 "/images/Horizontal_Axis_Camera/in_situ_absorption", ...
                 "/images/Vertical_Axis_Camera/in_situ_absorption"];

    % --- Validate camera index ---
    if options.cam < 1 || options.cam > numel(groupList)
        error('Invalid camera index: %d', options.cam);
    end

    files = dir(fullfile(options.folderPath, '*.h5'));
    nFiles = numel(files);
    if nFiles == 0
        error('No HDF5 files found in %s', options.folderPath);
    end

    % Determine image size
    testFile = fullfile(files(1).folder, files(1).name);
    atm_test = double(imrotate(h5read(testFile, append(groupList(options.cam), "/atoms")), options.angle, 'bilinear', 'crop'));
    [ny, nx] = size(atm_test);

    raw_scan_parameter_values = zeros(1, nFiles);
    raw_file_list             = string(zeros(1,nFiles)); % always string array
    
    if options.SAVE_TO_WORKSPACE
        fprintf('\n[INFO] Creating in-memory arrays of raw data...\n');
        full_od_imgs  = nan(ny, nx, nFiles, 'single');
        full_bkg_imgs = nan(ny, nx, nFiles, 'single');
    else
        % --- Create uniquely identified full OD image folder ---
        dataSource = makeDataSourceStruct(options.folderPath);
        runID = sprintf('%s_%s_Run%04d', ...
                        dataSource{1}.sequence, ...
                        strrep(dataSource{1}.date,'/','-'), ...
                        dataSource{1}.runs);
    
        % --- Determine parent folder for FullODImages ---
        if isfield(options, 'FullODImagesFoldersPath') && ...
           ~isempty(options.FullODImagesFoldersPath) && ...
           isfolder(options.FullODImagesFoldersPath)
            parentFolder = options.FullODImagesFoldersPath;
        else
            parentFolder = options.saveDirectory;
        end
        
        % --- Create uniquely identified full OD image folder ---
        fullODImageFolder = fullfile(parentFolder, ['FullODImages_' runID]);
        if ~exist(fullODImageFolder,'dir'), mkdir(fullODImageFolder); end
        fprintf('\n[INFO] Creating folder of full OD images on disk: %s\n', fullODImageFolder);

        % --- Save metadata for this run ---
        metadata.options   = options;
        metadata.timestamp = datetime;  % still record analysis time
        metadata.runID     = runID;     % traceable to experiment run
        metadata.imageSize = [ny, nx];
        metadata.fileList  = string(arrayfun(@(f) fullfile(f.folder, f.name), files, 'UniformOutput', false));
        save(fullfile(fullODImageFolder,'metadata.mat'),'metadata','-v7.3');
    
        full_od_imgs  = fullODImageFolder;
        full_bkg_imgs = fullODImageFolder;
    end

    % --- Prepare file names ---
    fullFileNames = string(arrayfun(@(f) fullfile(f.folder, f.name), files, 'UniformOutput', false));

    % --- Check for Parallel Computing Toolbox ---
    useParallel = license('test','Distrib_Computing_Toolbox') && ~options.SAVE_TO_WORKSPACE;

    if useParallel
        fprintf('\n[INFO] Parallel Computing Toolbox detected. Using parallelization for raw data processing...\n');

        % --- Preallocate scan parameters and file list for parallel mode ---
        raw_scan_parameter_values = nan(1, nFiles, 'single');
        raw_file_list             = fullFileNames;

        % --- Parallel loop without progress bar ---
        parfor k = 1:nFiles
            [od_img, bkg_img, val] = readAndComputeOD(fullFileNames(k), options, groupList, ny, nx);
            writeFullODImagesToDisk(fullODImageFolder, od_img, bkg_img, val, fullFileNames(k), k);
            raw_scan_parameter_values(k) = val;
        end
    else
        % --- Standard for-loop with your progress bar ---
        showPB = isfield(options,'showProgressBar') && options.showProgressBar;
        if showPB && options.SAVE_TO_WORKSPACE
            pb = Helper.ProgressBar();
            pb.run('Progress: ');
        end

        for k = 1:nFiles
            [od_img, bkg_img, val] = readAndComputeOD(fullFileNames(k), options, groupList, ny, nx);

            if options.SAVE_TO_WORKSPACE
                full_od_imgs(:,:,k)  = single(od_img);
                full_bkg_imgs(:,:,k) = single(bkg_img);
            else
                writeFullODImagesToDisk(fullODImageFolder, od_img, bkg_img, val, fullFileNames(k), k);
            end
            raw_scan_parameter_values(k) = val;
            raw_file_list(k)             = fullFileNames(k);

            if showPB && options.SAVE_TO_WORKSPACE
                progressPercent = round(k/nFiles*100);
                pb.run(progressPercent);
            end
        end

        if showPB && options.SAVE_TO_WORKSPACE
            pb.run(' Done!');
        end
    end
end

%% --- Local helper functions ---
function [od_img, bkg_img, val] = readAndComputeOD(fullFileName, options, groupList, ny, nx)
    try
        atm_img  = double(imrotate(h5read(fullFileName, append(groupList(options.cam), "/atoms")), options.angle, 'bilinear', 'crop'));
        bkg_img  = double(imrotate(h5read(fullFileName, append(groupList(options.cam), "/background")), options.angle, 'bilinear', 'crop'));
        dark_img = double(imrotate(h5read(fullFileName, append(groupList(options.cam), "/dark")), options.angle, 'bilinear', 'crop'));
        od_img   = Helper.calculateODImage(atm_img, bkg_img, dark_img, options.ImagingMode, options.PulseDuration);
    catch
        warning('\nMissing data in %s, storing NaNs.', fullFileName);
        od_img  = nan(ny, nx);
        bkg_img = nan(ny, nx);
    end

    % --- Extract scan parameter safely ---
    try
        info      = h5info(fullFileName, '/globals');
        attrNames = string({info.Attributes.Name});
        if ismember(options.scan_parameter, attrNames)
            val = h5readatt(fullFileName, '/globals', options.scan_parameter);
            if strcmp(options.scan_parameter,'ps_rot_mag_fin_pol_angle')
                val = 180 - val;
            end
        else
            val = NaN;
        end
    catch
        val = NaN;
    end
end

function writeFullODImagesToDisk(fullODImageFolder, od_img, bkg_img, scan_val, file_name, idx)
% Writes a single OD/BKG image + scan parameter to a MAT file
    matFilePath = fullfile(fullODImageFolder, sprintf('Image_%04d.mat', idx));
    OD  = single(od_img);
    BKG = single(bkg_img);
    Scan = single(scan_val);
    File = string(file_name);
    save(matFilePath, 'OD','BKG','Scan','File','-v7.3');
end

function dataSource = makeDataSourceStruct(folderPath)
    % Split by file separators (handles / or \)
    parts = regexp(folderPath, '[\\/]', 'split');

    % Remove empty parts caused by leading slashes
    parts = parts(~cellfun('isempty', parts));

    % Extract sequence, date, and run number
    % Now the indices are correct:
    % parts = {'DyLabNAS', 'Data', 'StructuralPhaseTransition', '2025', '08', '13', '0062'}
    sequence = parts{3};       % "StructuralPhaseTransition"
    year     = parts{4};       % "2025"
    month    = parts{5};       % "08"
    day      = parts{6};       % "13"
    runStr   = parts{7};       % "0062"

    % Build date string
    dateStr = sprintf('%s/%s/%s', year, month, day);

    % Convert run string to number
    runNum = str2double(runStr);

    % Construct struct inside a cell array
    dataSource = {
        struct('sequence', sequence, ...
               'date', dateStr, ...
               'runs', runNum)
    };
end