function Volume_Slicer()


% =====================================================================
% Volume_Slicer.m  -  Oblique Volume Slicer Suite (v2, 2026)
%
% A MATLAB GUI application for interactive inspection of 3D volumes
% (NRRD, NHDR, NIfTI) and batch export of 2D slice datasets (PNG) along
% standard or oblique planes.
%
% The GUI has been conceived and developed by Lionel Cervera Gontard 
% with the help of third-party libraries and vibe coding using Claude AI from Anthropic.
% ---------------------------------------------------------------------
% THIRD-PARTY CODE NOTICE
% ---------------------------------------------------------------------
% This file is self-contained: in order to avoid external dependencies,
% the NRRD/NHDR reader and its auxiliary helpers have been embedded as
% local functions at the bottom of this file. Full credit for those
% functions goes to their original author; they are redistributed here
% verbatim (modulo line-ending normalisation), under their original
% terms, and no algorithmic change has been introduced.
%
% Integrated third-party functions:
%
%   * nhdr_nrrd_read
%   * nrrd_getMatlabDataType
%   * nrrd_getSpaceDimensions
%     ... together with their internal sub-functions
%     (foundKeyword, extractStringList, extractNumbersWithout,
%      extract_spaceunits_list, extract_datafiles, extractGradient,
%      extractKeyValueString, readData, adjustEndian, sizeOf).
%
%   Original author : Gaetan Rensonnet
%   Original date   : July 2018
%   Original source : MATLAB File Exchange / NHDR-NRRD MATLAB toolbox
%
% The original reader module itself acknowledges prior contributions
% from Marc / K. Quintus (maReadNrrdHeader, maReadNrrdData,
% Harvard BWH), Jeff Mather (File Exchange NRRD reader) and
% "mdcacio" (File Exchange nrrdWriter); those attributions are
% preserved in the comment headers of the integrated functions below.
%
% All remaining code in this file (GUI, projection / rotation logic,
% intensity normalization, dataset export, imrotate_fill wrapper, etc.)
% is original work for the Oblique Volume Slicer Suite.
%
% Users who wish to cite, reuse or redistribute any of the integrated
% third-party functions should refer to their original publications /
% File Exchange entries and respect the licensing terms attached to
% them there.
% =====================================================================

    % Initialization
    gui = struct();
    gui.vol = [];
    gui.inputDir = ''; gui.outputDir = '';
    gui.currentFileName = '';
    gui.voxelSize = [NaN NaN NaN]; % [dx dy dz] display only
    gui.intStats = struct('min',NaN,'max',NaN,'p1',NaN,'p99',NaN);

    % --- GUI ---
    gui.fig = figure('Name', 'Oblique Volume Slicer Suite - v2 2026', 'NumberTitle', 'off', ...
        'Position', [50 50 1200 750], 'Color', [0.98 0.98 0.98], 'MenuBar', 'none');

    % =====================================================================
    % HELP MENU (NEW)
    % =====================================================================
    mHelp = uimenu(gui.fig, 'Label', 'Help');
    uimenu(mHelp, 'Label', 'Quick start', 'Callback', @showQuickStart);
    uimenu(mHelp, 'Label', 'Controls reference', 'Callback', @showControls);
    uimenu(mHelp, 'Label', 'About', 'Separator', 'on', 'Callback', @showAbout);

    % --- CONTROL PANEL ---
    gui.pCtrl = uipanel('Parent', gui.fig, 'Title', ' ⚙️ SETTINGS ', ...
        'Position', [0.01 0.02 0.22 0.96], 'BackgroundColor', 'white', 'FontWeight', 'bold');

    % Folders
    uicontrol('Parent', gui.pCtrl, 'Style', 'pushbutton', 'String', '📂 Input Folder', ...
        'Units', 'normalized', 'Position', [0.1 0.92 0.8 0.04], 'Callback', @selectIn);
    gui.lblIn = uicontrol('Parent', gui.pCtrl, 'Style', 'text', 'String', '...', ...
        'Units', 'normalized', 'Position', [0.1 0.89 0.8 0.02], 'FontSize', 8, 'BackgroundColor', 'white');

    uicontrol('Parent', gui.pCtrl, 'Style', 'pushbutton', 'String', '💾 Output Folder', ...
        'Units', 'normalized', 'Position', [0.1 0.82 0.8 0.04], 'Callback', @selectOut);
    gui.lblOut = uicontrol('Parent', gui.pCtrl, 'Style', 'text', 'String', '...', ...
        'Units', 'normalized', 'Position', [0.1 0.79 0.8 0.02], 'FontSize', 8, 'BackgroundColor', 'white');

    % File list
    gui.lstFiles = uicontrol('Parent', gui.pCtrl, 'Style', 'listbox', 'Units', 'normalized', ...
        'Position', [0.1 0.62 0.8 0.14], 'Callback', @loadFile);

    % Voxel size label
    gui.lblVoxel = uicontrol('Parent', gui.pCtrl, 'Style', 'text', ...
        'String', 'Voxel size: n/a', ...
        'Units', 'normalized', 'Position', [0.1 0.595 0.8 0.02], ...
        'FontSize', 9, 'BackgroundColor', 'white', 'HorizontalAlignment', 'left');

    % Volume slices label: Sagittal x Coronal x Axial
    gui.lblSlices = uicontrol('Parent', gui.pCtrl, 'Style', 'text', ...
        'String', 'Slices (S×C×A): n/a', ...
        'Units', 'normalized', 'Position', [0.1 0.575 0.8 0.02], ...
        'FontSize', 9, 'BackgroundColor', 'white', 'HorizontalAlignment', 'left');

    % Projection
    uicontrol('Parent', gui.pCtrl, 'Style', 'text', 'String', 'Projection:', 'Units', 'normalized', ...
        'Position', [0.1 0.55 0.8 0.02], 'BackgroundColor', 'white');
    gui.popProj = uicontrol('Parent', gui.pCtrl, 'Style', 'popupmenu', ...
        'String', {'Sagittal (X)', 'Coronal (Z)', 'Axial (Y)'}, ...
        'Units', 'normalized', 'Position', [0.1 0.52 0.8 0.03], 'Callback', @updateSliderRanges);

    % Rebin button
    gui.btnRebin = uicontrol('Parent', gui.pCtrl, 'Style', 'pushbutton', ...
        'String', '⬇ Rebin ×2', ...
        'Units', 'normalized', ...
        'Position', [0.1 0.49 0.8 0.035], ...
        'BackgroundColor', [0.85 0.85 0.85], ...
        'FontWeight', 'bold', ...
        'Enable', 'off', ...
        'Callback', @rebinVolume);

    % Hinge (axial only)
    gui.lblHinge = uicontrol('Parent', gui.pCtrl, 'Style', 'text', ...
        'String', 'Hinge axis (Axial):', 'Units', 'normalized', ...
        'Position', [0.1 0.47 0.8 0.02], 'BackgroundColor', 'white', 'Enable', 'off');
    gui.popHinge = uicontrol('Parent', gui.pCtrl, 'Style', 'popupmenu', ...
        'String', {'Hinge in Coronal Plane', 'Hinge in Sagittal Plane'}, ...
        'Units', 'normalized', 'Position', [0.1 0.44 0.8 0.03], ...
        'Enable', 'off', 'Callback', @updateUI);

    % Rotation axis position (+ NEW: No rotation)
    gui.lblAxisPos = uicontrol('Parent', gui.pCtrl, 'Style', 'text', ...
        'String', 'Rotation axis position:', ...
        'Units', 'normalized', 'Position', [0.1 0.41 0.8 0.02], ...
        'BackgroundColor', 'white');

    gui.popAxisPos = uicontrol('Parent', gui.pCtrl, 'Style', 'popupmenu', ...
        'String', {'No rotation', 'Slice center', 'Slice edge'}, ...
        'Units', 'normalized', 'Position', [0.1 0.38 0.8 0.03], ...
        'Value', 1, ...
        'Callback', @updateSliderRanges);

    % Angles (only main)
    uicontrol('Parent', gui.pCtrl, 'Style', 'text', 'String', 'Main angle:', 'Units', 'normalized', ...
        'Position', [0.1 0.35 0.8 0.02], 'BackgroundColor', 'white');
    gui.edAng = uicontrol('Parent', gui.pCtrl, 'Style', 'edit', 'String', '-21:5:20', ...
        'Units', 'normalized', 'Position', [0.1 0.32 0.8 0.03], 'Callback', @updateSliderRanges);

    % Export format (grayscale vs binary)
    uicontrol('Parent', gui.pCtrl, 'Style', 'text', 'String', 'Export output:', ...
        'Units', 'normalized', 'Position', [0.1 0.295 0.8 0.02], 'BackgroundColor', 'white', ...
        'HorizontalAlignment', 'left');
    gui.popOutType = uicontrol('Parent', gui.pCtrl, 'Style', 'popupmenu', ...
        'String', {'Grayscale (uint8 PNG)','Binary mask (0/255 PNG)'}, ...
        'Units', 'normalized', 'Position', [0.1 0.265 0.8 0.03], ...
        'Callback', @updateUI);

    % =====================================================================
    % NEW: Intensity normalization toggle (preview + grayscale export)
    % Checkbox ON: original behavior (p1–p99 stable normalization)
    % Checkbox OFF: preview autoscale, export uses global min/max
    % =====================================================================
    gui.chkNorm = uicontrol('Parent', gui.pCtrl, 'Style', 'checkbox', ...
        'String', 'Intensity normalization (p1–p99)', ...
        'Units', 'normalized', 'Position', [0.1 0.24 0.8 0.022], ...
        'Value', 1, ...
        'BackgroundColor', 'white', ...
        'Callback', @updateUI);

    % Output size (LAYOUT TWEAK: moved down so checkbox is not hidden)
    uicontrol('Parent', gui.pCtrl, 'Style', 'text', 'String', 'Output width (px):', ...
        'Units', 'normalized', 'Position', [0.1 0.21 0.35 0.02], 'BackgroundColor', 'white');
    gui.edW = uicontrol('Parent', gui.pCtrl, 'Style', 'edit', 'String', '256', ...
        'Units', 'normalized', 'Position', [0.1 0.172 0.35 0.03]);

    uicontrol('Parent', gui.pCtrl, 'Style', 'text', 'String', 'Output height (px):', ...
        'Units', 'normalized', 'Position', [0.55 0.21 0.35 0.02], 'BackgroundColor', 'white');
    gui.edH = uicontrol('Parent', gui.pCtrl, 'Style', 'edit', 'String', '256', ...
        'Units', 'normalized', 'Position', [0.55 0.172 0.35 0.03]);

    % --- GLOBAL SLICE RANGE (From / To / Step) ---
    uicontrol('Parent', gui.pCtrl, 'Style', 'text', ...
        'String', 'Export slice range (global):', ...
        'Units', 'normalized', 'Position', [0.1 0.145 0.8 0.02], ...
        'BackgroundColor', 'white', 'HorizontalAlignment', 'left');

    gui.edSliceFrom = uicontrol('Parent', gui.pCtrl, 'Style', 'edit', 'String', '1', ...
        'Units', 'normalized', 'Position', [0.1 0.115 0.28 0.03], ...
        'Callback', @clampGlobalRange);

    gui.edSliceTo = uicontrol('Parent', gui.pCtrl, 'Style', 'edit', 'String', '1', ...
        'Units', 'normalized', 'Position', [0.41 0.115 0.28 0.03], ...
        'Callback', @clampGlobalRange);

    gui.edSliceStep = uicontrol('Parent', gui.pCtrl, 'Style', 'edit', 'String', '1', ...
        'Units', 'normalized', 'Position', [0.72 0.115 0.18 0.03], ...
        'Callback', @clampGlobalRange);

    uicontrol('Parent', gui.pCtrl, 'Style', 'text', ...
        'String', 'From      To      Step', ...
        'Units', 'normalized', 'Position', [0.1 0.095 0.8 0.02], ...
        'BackgroundColor', 'white', 'HorizontalAlignment', 'left', 'FontSize', 8);

    gui.btnFullRange = uicontrol('Parent', gui.pCtrl, 'Style', 'pushbutton', ...
        'String', 'Full range', ...
        'Units', 'normalized', 'Position', [0.1 0.07 0.8 0.025], ...
        'Enable', 'off', 'Callback', @setFullGlobalRange);

    % Export button
    gui.btnProcess = uicontrol('Parent', gui.pCtrl, 'Style', 'pushbutton', 'String', '⚡ EXPORT DATASET', ...
        'Units', 'normalized', 'Position', [0.1 0.01 0.8 0.055], 'BackgroundColor', [0 0.4 0.5], ...
        'ForegroundColor', 'white', 'FontWeight', 'bold', 'Enable', 'off', 'Callback', @processDataset);

    % --- VIEW AREA ---
    gui.axStd = axes('Parent', gui.fig, 'Position', [0.28 0.55 0.32 0.38]);
    gui.sldStd = uicontrol('Parent', gui.fig, 'Style', 'slider', 'Units', 'normalized', ...
        'Position', [0.28 0.51 0.32 0.02], 'Callback', @updateUI);

    gui.axObl = axes('Parent', gui.fig, 'Position', [0.65 0.55 0.32 0.38]);
    gui.sldAng = uicontrol('Parent', gui.fig, 'Style', 'slider', 'Units', 'normalized', ...
        'Position', [0.65 0.51 0.32 0.02], 'Callback', @updateUI);

    gui.ax3D = axes('Parent', gui.fig, 'Position', [0.35 0.05 0.60 0.40]);

    % =====================================================================
    % FUNCTIONS
    % =====================================================================

    % ---------------- HELP menu callbacks (NEW) ----------------
    function showQuickStart(~, ~)
        msg = sprintf([ ...
            'Quick start:\n\n' ...
            '1) 📂 Input Folder: select a folder with .nrrd/.nhdr/.nii/.nii.gz\n' ...
            '2) 💾 Output Folder: choose where PNGs will be written\n' ...
            '3) Select a file from the list\n' ...
            '4) Choose Projection + (optional) rotation settings\n' ...
            '5) Choose Export output type (Grayscale or Binary mask)\n' ...
            '6) Intensity normalization:\n' ...
            '   - ON  : p1–p99 stable normalization (original behavior)\n' ...
            '   - OFF : preview autoscale (raw); export uses global min/max\n' ...
            '7) Set slice From/To/Step and click ⚡ EXPORT DATASET\n' ...
            ]);
        helpdlg(msg, 'Quick start');
    end

    function showControls(~, ~)
        msg = sprintf([ ...
            'Controls reference:\n\n' ...
            'Projection:\n' ...
            '  - Sagittal (X), Coronal (Z), Axial (Y)\n\n' ...
            'Rotation axis position:\n' ...
            '  - No rotation / Slice center / Slice edge\n\n' ...
            'Main angle:\n' ...
            '  - MATLAB range expression (e.g. -21:5:20)\n\n' ...
            'Export output:\n' ...
            '  - Grayscale (uint8 PNG)\n' ...
            '  - Binary mask (0/255 PNG)\n\n' ...
            'Intensity normalization:\n' ...
            '  - ON: p1–p99 scaling\n' ...
            '  - OFF: preview autoscale; export min/max\n' ...
            ]);
        helpdlg(msg, 'Controls reference');
    end

    function showAbout(~, ~)
        msg = sprintf([ ...
            'Oblique Volume Slicer Suite - v2 2026\n\n' ...
            'Added in this version:\n' ...
            '  - Help menu\n' ...
            '  - Intensity normalization checkbox\n' ...
            '  - Small layout tweak to avoid overlap\n' ...
            ]);
        msgbox(msg, 'About', 'help');
    end

    function selectIn(~, ~)
        d = uigetdir();
        if d ~= 0
            gui.inputDir = d;
            gui.lblIn.String = d;
            refreshFiles();
        end
    end

    function selectOut(~, ~)
        d = uigetdir();
        if d ~= 0
            gui.outputDir = d;
            gui.lblOut.String = d;
            checkReady();
        end
    end

    function refreshFiles()
        f = [ ...
            dir(fullfile(gui.inputDir, '*.nrrd')); ...
            dir(fullfile(gui.inputDir, '*.nhdr')); ...
            dir(fullfile(gui.inputDir, '*.nii'));  ...
            dir(fullfile(gui.inputDir, '*.nii.gz')) ...
        ];
        if ~isempty(f)
            gui.lstFiles.String = {f.name};
            gui.lstFiles.Value = 1;
        else
            gui.lstFiles.String = {'No files found'};
        end
    end

    function loadFile(~, ~)
        if isempty(gui.lstFiles.String) || startsWith(gui.lstFiles.String{1}, 'No')
            return;
        end

        gui.currentFileName = gui.lstFiles.String{gui.lstFiles.Value};
        fullFile = fullfile(gui.inputDir, gui.currentFileName);

        try
            [~,~,ext] = fileparts(gui.currentFileName);
            isNiiGz = endsWith(lower(gui.currentFileName), '.nii.gz');
            isNifti = false;

            gui.voxelSize = [NaN NaN NaN];
            set(gui.lblVoxel,  'String', 'Voxel size: n/a');
            set(gui.lblSlices, 'String', 'Slices (S×C×A): n/a');

            if strcmpi(ext, '.nrrd') || strcmpi(ext, '.nhdr')
                res = nhdr_nrrd_read(fullFile, true);
                vol = res.data;

                dx = NaN; dy = NaN; dz = NaN;
                try
                    if isfield(res,'meta')
                        m = res.meta;
                        if isfield(m,'spacings')
                            s = m.spacings;
                            if numel(s) >= 3, dx=s(1); dy=s(2); dz=s(3); end
                        elseif isfield(m,'space_directions')
                            sd = m.space_directions;
                            if isnumeric(sd) && all(size(sd) >= [3 3])
                                dx = norm(sd(:,1)); dy = norm(sd(:,2)); dz = norm(sd(:,3));
                            end
                        end
                    end
                catch
                end
                gui.voxelSize = [dx dy dz];

            elseif strcmpi(ext, '.nii')
                info = niftiinfo(fullFile);
                vol = niftiread(info);
                isNifti = true;
                if isfield(info,'PixelDimensions') && numel(info.PixelDimensions) >= 3
                    gui.voxelSize = double(info.PixelDimensions(1:3));
                end

            elseif isNiiGz
                tmpDir = tempname;
                mkdir(tmpDir);
                gunzip(fullFile, tmpDir);

                niiFiles = dir(fullfile(tmpDir, '*.nii'));
                if isempty(niiFiles)
                    error('Could not decompress .nii.gz (no .nii produced).');
                end

                niiPath = fullfile(tmpDir, niiFiles(1).name);
                info = niftiinfo(niiPath);
                vol = niftiread(info);
                isNifti = true;
                if isfield(info,'PixelDimensions') && numel(info.PixelDimensions) >= 3
                    gui.voxelSize = double(info.PixelDimensions(1:3));
                end

                delete(niiPath);
                rmdir(tmpDir);
            else
                error('Unsupported format: %s', gui.currentFileName);
            end

            if ndims(vol) > 3
                vol = vol(:,:,:,1);
            end

            if isNifti
                vol = permute(vol, [2 3 1]); % (Z,X,Y) -> (X,Y,Z)
            end

            if ~isfloat(vol)
                vol = single(vol);
            end

            gui.vol = vol;

            % Intensity stats (for stable preview/export scaling)
            try
                v = single(gui.vol(:));
                gui.intStats.min = min(v);
                gui.intStats.max = max(v);
                gui.intStats.p1  = prctile(v, 1);
                gui.intStats.p99 = prctile(v, 99);
            catch
                gui.intStats = struct('min',NaN,'max',NaN,'p1',NaN,'p99',NaN);
            end

            if all(isfinite(gui.voxelSize))
                set(gui.lblVoxel, 'String', sprintf('Voxel size: %.4g x %.4g x %.4g', ...
                    gui.voxelSize(1), gui.voxelSize(2), gui.voxelSize(3)));
            else
                set(gui.lblVoxel, 'String', 'Voxel size: n/a');
            end

            nx = size(gui.vol,1);
            ny = size(gui.vol,2);
            nz = size(gui.vol,3);
            set(gui.lblSlices, 'String', sprintf('Slices (S×C×A): %d × %d × %d', nx, nz, ny));

            set(gui.btnRebin, 'Enable', 'on');
            set(gui.btnFullRange, 'Enable', 'on');

            if isempty(strtrim(gui.edSliceFrom.String)) || isempty(strtrim(gui.edSliceTo.String)) || ...
               strcmp(strtrim(gui.edSliceFrom.String),'1') && strcmp(strtrim(gui.edSliceTo.String),'1')
                setFullGlobalRange();
            end
            clampGlobalRange();

            updateSliderRanges();
            checkReady();

        catch ME
            errordlg(ME.message, 'Volume load error');
        end
    end

    % ---------------- GLOBAL RANGE helpers ----------------
    function maxGlobal = getMaxGlobalSlices()
        if isempty(gui.vol)
            maxGlobal = 1;
            return;
        end
        s = size(gui.vol);
        if numel(s) < 3
            maxGlobal = max(s);
        else
            maxGlobal = max(s(1:3));
        end
    end

    function [iFrom, iTo, iStep] = getGlobalRange()
        iFrom = round(str2double(gui.edSliceFrom.String));
        iTo   = round(str2double(gui.edSliceTo.String));
        iStep = round(str2double(gui.edSliceStep.String));

        if isnan(iFrom), iFrom = 1; end
        if isnan(iTo),   iTo = 1; end
        if isnan(iStep) || iStep < 1, iStep = 1; end

        if iFrom > iTo
            tmp = iFrom; iFrom = iTo; iTo = tmp;
        end

        iFrom = max(1, iFrom);
        iTo   = max(1, iTo);
    end

    function clampGlobalRange(~, ~)
        [iFrom, iTo, iStep] = getGlobalRange();
        set(gui.edSliceFrom, 'String', num2str(iFrom));
        set(gui.edSliceTo,   'String', num2str(iTo));
        set(gui.edSliceStep, 'String', num2str(max(1, iStep)));
    end

    function setFullGlobalRange(~, ~)
        if isempty(gui.vol), return; end
        maxG = getMaxGlobalSlices();
        set(gui.edSliceFrom, 'String', '1');
        set(gui.edSliceTo,   'String', num2str(maxG));
        set(gui.edSliceStep, 'String', '1');
    end

    % ---------------- Rebin ----------------
    function rebinVolume(~, ~)
        if isempty(gui.vol)
            return;
        end

        set(gui.btnRebin, 'String', 'Processing...', 'Enable', 'off');
        drawnow;

        try
            vol = gui.vol;
            [nx, ny, nz] = size(vol);

            nx2 = max(1, floor(nx/2));
            ny2 = max(1, floor(ny/2));
            nz2 = max(1, floor(nz/2));

            vol2 = zeros(nx2, ny2, nz2, 'like', vol);

            for x = 1:nx2
                for y = 1:ny2
                    for z = 1:nz2
                        bx = (2*x-1):min(2*x, nx);
                        by = (2*y-1):min(2*y, ny);
                        bz = (2*z-1):min(2*z, nz);

                        block = vol(bx, by, bz);
                        vol2(x,y,z) = mean(block(:), 'omitnan');
                    end
                end
            end

            gui.vol = vol2;

            if all(isfinite(gui.voxelSize))
                gui.voxelSize = gui.voxelSize * 2;
                set(gui.lblVoxel, 'String', sprintf('Voxel size: %.4g x %.4g x %.4g', ...
                    gui.voxelSize(1), gui.voxelSize(2), gui.voxelSize(3)));
            end

            nx = size(gui.vol,1);
            ny = size(gui.vol,2);
            nz = size(gui.vol,3);
            set(gui.lblSlices, 'String', sprintf('Slices (S×C×A): %d × %d × %d', nx, nz, ny));

            clampGlobalRange();

            updateSliderRanges();
            updateUI();

        catch ME
            errordlg(ME.message, 'Rebin error');
        end

        set(gui.btnRebin, 'String', '⬇ Rebin ×2', 'Enable', 'on');
        drawnow;
    end

    % ---------------- Sliders / UI ----------------
    function updateSliderRanges(~, ~)
        if isempty(gui.vol), return; end
        [nx, ny, nz] = size(gui.vol);
        mode = gui.popProj.Value;

        % Hinge only for axial (and only meaningful if rotation is ON)
        axisOpt = gui.popAxisPos.Value; % 1=No rot, 2=center, 3=edge
        rotOn = (axisOpt ~= 1);

        st = 'off';
        if (mode == 3) && rotOn, st = 'on'; end
        set(gui.lblHinge, 'Enable', st);
        set(gui.popHinge, 'Enable', st);

        % Disable main angle controls if rotation is off
        set(gui.edAng,  'Enable', ternary(rotOn,'on','off'));
        set(gui.sldAng, 'Enable', ternary(rotOn,'on','off'));

        dims = [nx, nz, ny];
        maxS = dims(mode);
        set(gui.sldStd, 'Min', 1, 'Max', maxS, 'Value', round(maxS/2));

        setupSnapSlider(gui.sldAng, gui.edAng.String);

        updateUI();
    end

    function out = ternary(cond, a, b)
        if cond, out = a; else, out = b; end
    end

    function setupSnapSlider(sld, strRange)
        try
            r = eval(strRange);
            if length(r) > 1
                span = max(r)-min(r);
                s = r(2)-r(1);
                set(sld, 'Min', min(r), 'Max', max(r), 'Value', r(1), ...
                    'SliderStep', [abs(s/span), abs(s/span)]);
            else
                set(sld, 'Min', r(1)-0.1, 'Max', r(1)+0.1, 'Value', r(1));
            end
        catch
            set(sld, 'Min', -1, 'Max', 1, 'Value', 0);
        end
    end

    function imgD = normalizeForDisplay(img)
        % Stable display mapping so background stays black across angles
        img = single(img);
        if isfield(gui,'popOutType') && gui.popOutType.Value == 2
            % binary preview (assumes mask is 0/1 or 0/255)
            if ~isfinite(gui.intStats.max) || gui.intStats.max == 0
                thr = 0.5;
            else
                thr = 0.5 * gui.intStats.max;
            end
            imgD = single(img > thr);
            return;
        end

        % NEW: allow disabling normalization for preview
        normOn = true;
        if isfield(gui,'chkNorm') && ishghandle(gui.chkNorm)
            normOn = logical(gui.chkNorm.Value);
        end
        if ~normOn
            imgD = img; % raw (imagesc autoscale handled in updateUI)
            return;
        end

        % Grayscale preview: use volume-global percentiles
        lo = gui.intStats.p1;
        hi = gui.intStats.p99;
        if ~isfinite(lo) || ~isfinite(hi) || hi <= lo
            lo = gui.intStats.min;
            hi = gui.intStats.max;
        end
        % clip to suppress bicubic overshoot / negatives
        img = min(max(img, gui.intStats.min), gui.intStats.max);
        imgD = (img - lo) ./ max(eps, (hi - lo));
        imgD = min(max(imgD, 0), 1);
    end

    function updateUI(~, ~)
        if isempty(gui.vol), return; end

        axisOpt  = gui.popAxisPos.Value; % 1=No rot, 2=center, 3=edge
        angM = 0;
        if axisOpt ~= 1
            angM = getSnappedVal(gui.sldAng, gui.edAng.String);
        end

        sliceIdx = round(gui.sldStd.Value);
        mode     = gui.popProj.Value;
        hType    = gui.popHinge.Value;

        imgStd = extractSlice(gui.vol, sliceIdx, mode);

        % ===========================
        % ROTATION LOGIC (UNCHANGED)
        % ===========================
        if axisOpt == 1
            vRot = gui.vol;
        else
            axisPos = axisOpt - 1; % 2->1(center), 3->2(edge)
            if isfield(gui,'popOutType') && gui.popOutType.Value == 2
                interp = 'nearest';
            else
                interp = 'bilinear';
            end
            vRot = rotateVolumeSingle(gui.vol, angM, mode, hType, axisPos, interp, 0);
        end
        % ===========================

        imgObl = extractSlice(vRot, sliceIdx, mode);

        % NEW: if normalization is OFF, preview uses imagesc autoscale
        normOn = true;
        if isfield(gui,'chkNorm') && ishghandle(gui.chkNorm)
            normOn = logical(gui.chkNorm.Value);
        end
        isBinary = (isfield(gui,'popOutType') && gui.popOutType.Value == 2);

        axes(gui.axStd);
        imgStdD = normalizeForDisplay(imgStd);
        if isBinary || normOn
            imagesc(imgStdD,[0 1]); colormap gray; axis image;
        else
            imagesc(imgStdD); colormap gray; axis image;
        end
        set(gui.axStd, 'XTick', [], 'YTick', []);

        maxS = round(get(gui.sldStd,'Max'));
        switch mode
            case 1, viewName = 'Standard view (Sagittal)';
            case 2, viewName = 'Standard view (Coronal)';
            case 3, viewName = 'Standard view (Axial)';
        end
        title(sprintf('%s | Slice %d / %d', viewName, sliceIdx, maxS));

        axes(gui.axObl);
        imgOblD = normalizeForDisplay(imgObl);
        if isBinary || normOn
            imagesc(imgOblD,[0 1]); colormap gray; axis image;
        else
            imagesc(imgOblD); colormap gray; axis image;
        end
        set(gui.axObl, 'XTick', [], 'YTick', []);

        if axisOpt == 1
            title(gui.axObl, 'Oblique: (rotation OFF)');
        else
            title(gui.axObl, sprintf('Oblique: %+d°', round(angM)));
        end

        draw3DPreview(sliceIdx, mode, angM, hType, axisOpt);
    end

    function draw3DPreview(sliceIdx, mode, angM, hType, axisOpt)
        % axisOpt: 1=No rot, 2=center, 3=edge

        axes(gui.ax3D);
        cla; hold on;
        [nx, ny, nz] = size(gui.vol);

        set(gui.ax3D,'DataAspectRatio',[1 1 1]);
        xlim([0 nx]); ylim([0 ny]); zlim([0 nz]);
        grid on; box on; view(135, 30);
        set(gui.ax3D, 'XTick', [], 'YTick', [], 'ZTick', []);

        vertices = [0 0 0; nx 0 0; nx ny 0; 0 ny 0; 0 0 nz; nx 0 nz; nx ny nz; 0 ny nz];
        faces = [1 2 6 5; 2 3 7 6; 3 4 8 7; 4 1 5 8; 1 2 3 4; 5 6 7 8];
        patch('Vertices', vertices, 'Faces', faces, 'FaceColor', 'none', ...
              'EdgeColor', [0.6 0.6 0.6], 'LineStyle', '--');

        nGrid = 10;
        switch mode
            case 1
                [Y, Z] = meshgrid(linspace(0, ny, nGrid), linspace(0, nz, nGrid));
                X = sliceIdx * ones(size(Y));
                projLabel = 'SAGITTAL';
            case 2
                [X, Y] = meshgrid(linspace(0, nx, nGrid), linspace(0, ny, nGrid));
                Z = sliceIdx * ones(size(X));
                projLabel = 'CORONAL';
            case 3
                [X, Z] = meshgrid(linspace(0, nx, nGrid), linspace(0, nz, nGrid));
                Y = sliceIdx * ones(size(X));
                projLabel = 'AXIAL';
        end

        if axisOpt == 1
            surf(X, Y, Z, 'FaceAlpha', 0.3, 'EdgeAlpha', 0.2);
            text(mean(X(:)), mean(Y(:)), mean(Z(:)), [projLabel ' (OFF)'], ...
                'HorizontalAlignment','center', 'FontWeight', 'bold', ...
                'BackgroundColor', 'yellow', 'EdgeColor', 'black');
            hold off;
            return;
        end

        axisPos = axisOpt - 1; % 1=center, 2=edge

        radM = -deg2rad(angM);

        rotX = @(t) [1 0 0; 0 cos(t) -sin(t); 0 sin(t) cos(t)];
        rotZ = @(t) [cos(t) -sin(t) 0; sin(t) cos(t) 0; 0 0 1];

        R = eye(3);
        switch mode
            case 1
                R = rotZ(radM);
            case 2
                R = rotX(radM);
            case 3
                if hType==1
                    R = rotX(radM);
                else
                    R = rotZ(radM);
                end
        end

        switch mode
            case 1
                dirWorld = [0 0 1];
                if axisPos == 1, p0 = [sliceIdx, ny/2, nz/2]; else, p0 = [sliceIdx, 0, nz/2]; end
            case 2
                dirWorld = [1 0 0];
                if axisPos == 1, p0 = [nx/2, ny/2, sliceIdx]; else, p0 = [nx/2, 0, sliceIdx]; end
            case 3
                if hType == 1
                    dirWorld = [1 0 0];
                    if axisPos == 1, p0 = [nx/2, sliceIdx, nz/2]; else, p0 = [nx/2, sliceIdx, 0]; end
                else
                    dirWorld = [0 0 1];
                    if axisPos == 1, p0 = [nx/2, sliceIdx, nz/2]; else, p0 = [0, sliceIdx, nz/2]; end
                end
        end

        pts = [X(:)'; Y(:)'; Z(:)'];
        pts_rel = pts - p0';
        pts_rot = R * pts_rel;
        pts_world = pts_rot + p0';

        Xf = reshape(pts_world(1,:), size(X));
        Yf = reshape(pts_world(2,:), size(Y));
        Zf = reshape(pts_world(3,:), size(Z));

        surf(Xf, Yf, Zf, 'FaceColor', 'cyan', 'FaceAlpha', 0.6, ...
            'EdgeColor', 'blue', 'EdgeAlpha', 0.3);

        text(mean(Xf(:)), mean(Yf(:)), mean(Zf(:)), projLabel, ...
            'HorizontalAlignment','center', 'FontWeight', 'bold', ...
            'BackgroundColor', 'yellow', 'EdgeColor', 'black');

        L = 1.2 * max([nx, ny, nz]);
        dirUnit = dirWorld / max(eps, norm(dirWorld));
        p1g = p0 - 0.5 * L * dirUnit;
        p2g = p0 + 0.5 * L * dirUnit;

        plot3([p1g(1) p2g(1)], [p1g(2) p2g(2)], [p1g(3) p2g(3)], 'r-', 'LineWidth', 2);
        pmid = (p1g + p2g) / 2;
        plot3(pmid(1), pmid(2), pmid(3), 'ro', 'MarkerFaceColor', 'r');

        hold off;
    end

    function val = getSnappedVal(sld, strRange)
        try
            r = eval(strRange);
            [~, idx] = min(abs(r - sld.Value));
            val = r(idx);
            set(sld, 'Value', val);
        catch
            val = 0;
        end
    end

    function slice = extractSlice(vol, idx, mode)
        [nx, ny, nz] = size(vol);
        try
            switch mode
                case 1
                    idx = min(max(1,idx),nx);
                    slice = squeeze(vol(idx, :, :));
                case 2
                    idx = min(max(1,idx),nz);
                    slice = squeeze(vol(:, :, idx))';
                case 3
                    idx = min(max(1,idx),ny);
                    slice = squeeze(vol(:, idx, :))';
            end
        catch
            slice = zeros(10,10);
        end
    end

    % ---------------------------------------------------------------------
    % REAL volume rotation
    % Single-angle version (fine adjustment removed)
    % axisPos: 1=center, 2=edge
    % ---------------------------------------------------------------------
    function vRot = rotateVolumeSingle(vol, angM, mode, hType, axisPos, interpMethod, fillVal)
        if nargin < 5
            axisPos = 1;
        end
        if nargin < 6 || isempty(interpMethod)
            interpMethod = 'bilinear';
        end
        if nargin < 7 || isempty(fillVal)
            fillVal = 0;
        end

        if angM == 0
            vRot = vol;
            return;
        end

        [nx, ny, nz] = size(vol);

        switch mode
            case 1 % Sagittal (X)
                if axisPos == 1
                    vRot = imrotate_fill(vol, angM, interpMethod, 'crop', fillVal);
                else
                    vP = cat(2, zeros(nx, ny, nz, 'like', vol), vol);
                    vR = imrotate_fill(vP, angM, interpMethod, 'crop', fillVal);
                    vRot = vR(:, ny+1:end, :);
                end

            case 2 % Coronal (Z)
                if axisPos == 1
                    vPerm = permute(vol, [2 3 1]);      % (y,z,x)
                    vR    = imrotate_fill(vPerm, angM, interpMethod, 'crop', fillVal);
                    vRot  = ipermute(vR, [2 3 1]);
                else
                    vPerm = permute(vol, [2 3 1]);      % (y,z,x)
                    vP    = cat(1, zeros(ny, nz, nx, 'like', vol), vPerm);
                    vR    = imrotate_fill(vP, angM, interpMethod, 'crop', fillVal);
                    vRotPerm = vR(ny+1:end, :, :);
                    vRot  = ipermute(vRotPerm, [2 3 1]);
                end

            case 3 % Axial (Y)
                if hType == 1
                    if axisPos == 1
                        vRotF = zeros(nx, ny, nz, 'like', vol);
                        for i = 1:nx
                            vRotF(i,:,:) = imrotate_fill(squeeze(vol(i,:,:)), angM, interpMethod, 'crop', fillVal);
                        end
                        vRot = vRotF;
                    else
                        vP = cat(3, zeros(nx, ny, nz, 'like', vol), vol);
                        vRotF = zeros(nx, ny, 2*nz, 'like', vol);
                        for i = 1:nx
                            vRotF(i,:,:) = imrotate_fill(squeeze(vP(i,:,:)), angM, interpMethod, 'crop', fillVal);
                        end
                        vRot = vRotF(:,:,nz+1:end);
                    end
                else
                    if axisPos == 1
                        vRotF = zeros(nx, ny, nz, 'like', vol);
                        for i = 1:nz
                            vRotF(:,:,i) = imrotate_fill(squeeze(vol(:,:,i)), angM, interpMethod, 'crop', fillVal);
                        end
                        vRot = vRotF;
                    else
                        vP = cat(1, zeros(nx, ny, nz, 'like', vol), vol);
                        vRotF = zeros(2*nx, ny, nz, 'like', vol);
                        for i = 1:nz
                            vRotF(:,:,i) = imrotate_fill(squeeze(vP(:,:,i)), angM, interpMethod, 'crop', fillVal);
                        end
                        vRot = vRotF(nx+1:end,:,:);
                    end
                end
        end
    end

    function checkReady()
        if ~isempty(gui.vol) && ~isempty(gui.outputDir)
            gui.btnProcess.Enable = 'on';
        end
    end

    function processDataset(~, ~)
    try
        if isempty(gui.vol)
            errordlg('No volume loaded.', 'Export error');
            return;
        end
        if isempty(gui.outputDir)
            errordlg('Please select an output folder.', 'Export error');
            return;
        end

        axisOpt = gui.popAxisPos.Value; % 1=No rot, 2=center, 3=edge

        if axisOpt == 1
            angsM = 0;
        else
            angsM = eval(gui.edAng.String);
        end

        mIdx  = gui.popProj.Value;
        hIdx  = gui.popHinge.Value;

        outW = str2double(gui.edW.String);
        outH = str2double(gui.edH.String);
        if isnan(outW) || outW <= 0, outW = []; end
        if isnan(outH) || outH <= 0, outH = []; end

        % --- folder tags ---
        projName = gui.popProj.String{mIdx};
        projTag  = regexprep(projName, '[^A-Za-z0-9]+', '');

        if axisOpt == 1
            axisTag = 'NoRotation';
            axisPos = 1; % dummy
        elseif axisOpt == 2
            axisTag = 'AxisCenter';
            axisPos = 1;
        else
            axisTag = 'AxisEdge';
            axisPos = 2;
        end

        if mIdx == 3 && axisOpt ~= 1
            if hIdx == 1
                hingeTag = 'HingeCoronal';
            else
                hingeTag = 'HingeSagittal';
            end
        else
            hingeTag = 'HingeNA';
        end

        [~, nameOnly] = fileparts(gui.currentFileName);
        fullP = fullfile(gui.outputDir, ...
            sprintf('%s_%s_%s_%s_Dataset', nameOnly, projTag, axisTag, hingeTag));

        if ~exist(fullP, 'dir')
            mkdir(fullP);
        end

        [gFrom, gTo, gStep] = getGlobalRange();
        gStep = max(1, gStep);

        total = length(angsM);
        if total == 0, total = 1; end

        h = waitbar(0, 'Exporting...');
        count = 0;

        for aM = angsM
            count = count + 1;
            waitbar(count/max(1,total), h);

            % ===========================
            % ROTATION LOGIC (UNCHANGED)
            % ===========================
            if axisOpt == 1
                vR = gui.vol;
            else
                if isfield(gui,'popOutType') && gui.popOutType.Value == 2
                interp = 'nearest';
            else
                interp = 'bilinear';
            end
            vR = rotateVolumeSingle(gui.vol, aM, mIdx, hIdx, axisPos, interp, 0);
            end
            % ===========================

            nS = size(vR, mIdx);
            if mIdx==2
                nS = size(vR,3);
            elseif mIdx==3
                nS = size(vR,2);
            end

            iFrom = max(1, min(nS, gFrom));
            iTo   = max(1, min(nS, gTo));
            if iFrom > iTo
                tmp = iFrom; iFrom = iTo; iTo = tmp;
            end

            for i = iFrom:gStep:iTo
                img = extractSlice(vR, i, mIdx);
                if isempty(img), continue; end

                outType = 1;
                if isfield(gui,'popOutType')
                    outType = gui.popOutType.Value; % 1=grayscale, 2=binary
                end

                if outType == 2
                    % --- Binary mask export (stable) ---
                    img = single(img);

                    % Heuristic threshold: assumes masks are 0/1 or 0/255. Use half max.
                    if isfinite(gui.intStats.max) && gui.intStats.max ~= 0
                        thr = 0.5 * gui.intStats.max;
                    else
                        thr = 0.5;
                    end
                    imgB = img > thr;

                    if ~isempty(outW) && ~isempty(outH)
                        imgB = imresize(imgB, [outH outW], 'nearest', 'Antialiasing', false);
                    end

                    imgU8 = uint8(255 * imgB);

                else
                    % --- Grayscale intensity export (stable, volume-global scaling) ---
                    img = single(img);

                    % NEW: checkbox controls export scaling
                    normOn = true;
                    if isfield(gui,'chkNorm') && ishghandle(gui.chkNorm)
                        normOn = logical(gui.chkNorm.Value);
                    end

                    if normOn
                        % ON: original p1/p99 scaling
                        lo = gui.intStats.p1;
                        hi = gui.intStats.p99;
                        if ~isfinite(lo) || ~isfinite(hi) || hi <= lo
                            lo = gui.intStats.min;
                            hi = gui.intStats.max;
                        end

                        % clip to suppress bicubic overshoot / negatives
                        img = min(max(img, gui.intStats.min), gui.intStats.max);

                        imgN = (img - lo) ./ max(eps, (hi - lo));
                        imgN = min(max(imgN, 0), 1);
                    else
                        % OFF: global min/max scaling (no percentile normalization/clipping)
                        lo = gui.intStats.min;
                        hi = gui.intStats.max;
                        if ~isfinite(lo) || ~isfinite(hi) || hi <= lo
                            lo = min(img(:));
                            hi = max(img(:));
                            if ~isfinite(lo) || ~isfinite(hi) || hi <= lo
                                lo = 0; hi = 1;
                            end
                        end

                        imgN = (img - lo) ./ max(eps, (hi - lo));
                        imgN = min(max(imgN, 0), 1);
                    end

                    if ~isempty(outW) && ~isempty(outH)
                        imgN = imresize(imgN, [outH outW], 'bilinear');
                    end

                    imgU8 = uint8(255 * imgN);
                end

                % File name without fine angle
                imwrite(repmat(imgU8,[1 1 3]), ...
                    fullfile(fullP, sprintf('Slice_%03d_P%+d.png', ...
                    i, round(aM))));
            end
        end

        close(h);
        msgbox('Dataset complete.');
    catch ME
        errordlg(ME.message, 'Export error');
    end
end

end

function B = imrotate_fill(A, ang, method, bbox, fillVal)
%IMROTATE_FILL Backward-compatible imrotate with configurable fill value.
% Some MATLAB versions do not support the 'FillValues' name-value pair in imrotate.
% This wrapper uses classic imrotate syntax and, when fillVal ~= 0, applies the fill
% by rotating a binary mask and assigning fillVal outside the original support.
%
% Supports 2-D arrays. If A is 3-D, rotates each A(:,:,k) plane independently.

    if nargin < 5 || isempty(fillVal)
        fillVal = 0;
    end
    if nargin < 4 || isempty(bbox)
        bbox = 'crop';
    end
    if nargin < 3 || isempty(method)
        method = 'nearest';
    end

    if ndims(A) > 2
        B = zeros(size(A), 'like', A);
        for k = 1:size(A,3)
            B(:,:,k) = imrotate_fill(A(:,:,k), ang, method, bbox, fillVal);
        end
        return;
    end

    % Classic, widely-supported imrotate signature
    B = imrotate(A, ang, method, bbox);

    % Apply non-zero fill if requested
    if fillVal ~= 0
        m = imrotate(ones(size(A), 'like', A), ang, 'nearest', bbox);
        B(m == 0) = fillVal;
    end
end

% =====================================================================
% Integrated external helpers (originally distributed as separate files)
% - nhdr_nrrd_read
% - nrrd_getMatlabDataType
% - nrrd_getSpaceDimensions
% and their internal sub-functions.
% Source: NHDR/NRRD Matlab toolbox by Gaetan Rensonnet (July 2018).
% =====================================================================

% Read NHDR/NRRD files 
% 
% headerInfo = nhdr_nrrd_read(nrrdFileName, bReadData)
%
% Inputs:
% * nrrdFileName (char): path to NHDR/NRRD file, either a detached header,
% a detached header pointing to detached data files or a NRRD standalone
% file with header and data included
% * bReadData (bool): set to true to read the data and store it in
% headerInfo.data, set to false to just import the header without the data
%
% Outputs:
% * headerInfo (struct): contains all the fields specified in the file
% header. If data was read it is contained in the 'data' field. That
% structure can be fed to the nhdr_nrrd_write module as is and produce
% valid NHDR/NRRD files.
% 
% Format definition:
% http://teem.sourceforge.net/nrrd/format.html
%
% A few supported NRRD features:
% - detached headers with all variants of 'data file:'
% - raw, txt/text/ascii, hex, gz/gzip encodings
% - definition of space and orientation
% - handling of diffusion-weighted MRI data with '<key>:=<value>' lines 
%   'modality:=DWMRI', 'DWMRI_b-value:=' and 'DWMRI_gradient_0000:=', 
%   'DWMRI_gradient_0001:=', 'DWMRI_gradient_0002:=', etc.
%   (see https://www.na-mic.org/wiki/NAMIC_Wiki:DTI:Nrrd_format)
% 
%
% Other features:
% - exits cleanly upon error (no file should be left open)
% - informative error messages whenever possible
% - warnings printed to console 
% - NHDR/NRRD writer fills out missing fields by inspecting the data
%   whenever possible
% 
%
% Unsupported features:
%
% In general,any header field that we were unable to parse is reported
% in a message printed to the console. Specific examples of
% unsupported features include:
%
% - reading data along more than 1 dimension or along a dimension other
%   than the slowest (last) axis specified by the optional <subdim>, as in
%   'data file: <format> <min> <max> <step> [<subdim>]' or 'data file: LIST
%   [<subdim>]'
% - storing all the comments found in headers
% - hex encoding
% - bz2/bzip2 encoding
% - byte skip; can only accept -1 with raw encoding, 0 for all other
%   encodings
% - line skip: can only accept 0
% - checking that field specifications of the form '<field>: <desc>' 
%   appear no more than once in the NRRD header (unlike '<key>:=<value>'
%   lines)
%
% Date: July 2018
% Author: Gaetan Rensonnet
%
%
%
% Notes :
%
% This is intended as a more comprehensive and accurate implementation of
% the NRRD format specification than most of the Matlab scripts that have
% been proposed so far. We try to fail graciously when an unsupported
% feature of the NRRD format is encountered. One of our main contributions
% is to propose a writer module which is compatible with the read module,
% in that the output of one can be given as an argument to the other to
% read or produce equivalent NHDR/NRRD files. This is still a version with
% much room for improvement.
%
% The following contributions inspired parts of our Matlab read/write
% modules:
% 
% 1. The body of the writer module was pretty much written from scratch but
% the general structure of the reader's main body is based on the Matlab
% functions maReadNrrdHeader.m and maReadNrrdData.m by marc@bwh.harvard.edu
% and kquintus@bwh.harvard.edu (unpublished). Many additions were made and
% a few bugs fixed.
%
% 2. Jeff Mather's NRRD reader
% (http://nl.mathworks.com/matlabcentral/fileexchange/34653-nrrd-format-file-reader)
% and
% http://jeffmatherphotography.com/dispatches/2012/02/writing-a-file-reader-in-matlab/)
% provided the auxiliary functions:
% - adjustEndian
% - getDatatype, which we renamed nrrd_getMatlabDataType and now throws
% a gracious error if it encounters 'block'-type data,
% - readData was adapted to include a cross-platform fix to delete
% temporary files when using gzip encoding. David Feng's fix used a
% Windows-specific command to delete temporary files
% (https://nl.mathworks.com/matlabcentral/fileexchange/50830-nrrd-format-file-reader).
%
% 3. mdcacio's nrrdWriter
% (https://nl.mathworks.com/matlabcentral/fileexchange/48621-nrrdwriter-filename--matrix--pixelspacing--origin--encoding-)
% provided the auxiliary functions:
% - writeData(): we got rid of the 'unexpected end of input stream when
% attempting to gunzip the file' error when using gzip encoding, which we
% later found had been fixed by Quan Chen independenlty and in a similar
% manner.
% - setDatatype(), which we renamed get_nrrd_datatype() and is the
% reciprocal of getDatatype 
%
function headerInfo = nhdr_nrrd_read(nrrdFileName, bReadData)
[mainFpath,mainFileName,mainFext] = fileparts( nrrdFileName );
headerInfo = struct();
% default header:
headerInfo.content = mainFileName;        % default value which can be overriden if content field is set
headerInfo.data = [];
fidr = fopen( nrrdFileName, 'r' );
if (fidr == -1)
    error('nhdr_nrrd_read ABORT: %s does not exist.\n', nrrdFileName );
end
try 
    if ~(strcmpi(mainFext, '.nhdr') || strcmpi(mainFext, '.nrrd'))
        warning('%s looks like a %s file, not a nhdr or nrrd file.\n', nrrdFileName, mainFext );
    end
    
    % Magic line
    cs  = fgetl(fidr);
    assert(numel(cs)>=4 && isequal(cs(1:4), 'NRRD'), 'Bad signature. First line should be magic line of type "NRRD000X" with X an integer between 1 and 5.'); % FIXME should throw an error if a bad integer is provided
    
    % Always optional: defining orientation
    define_orientation = 0;         % internal-use variable
    
    % Parse header
    while ( ~feof(fidr) )
        
        cs = fgetl( fidr );         % content string
        
        if isempty(cs)
            % End of header
            break;
        end
                
        if ( foundKeyword( 'CONTENT:', cs ) )
            
            headerInfo.content = strtrim( cs( length('CONTENT:')+1:end ) );
            
        elseif ( foundKeyword('TYPE:', cs ) )
            
            headerInfo.type = strtrim( cs( length('TYPE:')+1:end ) );
            
        elseif ( foundKeyword('ENDIAN:', cs ) )
            
            headerInfo.endian = strtrim( cs( length('ENDIAN:')+1:end ) );
            
        elseif ( foundKeyword('ENCODING:', cs ) )
            
            headerInfo.encoding = strtrim( cs( length('ENCODING:')+1:end ) );
            
        elseif ( foundKeyword('DIMENSION:', cs ) )
            
            headerInfo.dimension = sscanf( cs( length('DIMENSION:')+1:end ), '%i' );
            
        elseif ( foundKeyword('SIZES:', cs ) )
            
            iSizes = sscanf( cs(length('SIZES:')+1:end), '%i' );
            headerInfo.sizes = iSizes(:)';      % store as row vector
            
        elseif ( foundKeyword('KINDS:', cs ) )
            
            headerInfo.kinds = extractStringList( cs(length('KINDS:')+1:end) ); % bug fixed with extractStringList where 2 entries are present
            % FIXME: check that axis sizes match each kind according to standard
        elseif ( foundKeyword('SPACINGS:', cs ) ) 
                ind=find(cs==' ');
                spacings(1)=str2num(cs(ind(1):ind(2)));
                spacings(2)=str2num(cs(ind(2):ind(3)));
                spacings(3)=str2num(cs(ind(3):end));
                headerInfo.spacings=spacings;
        elseif ( foundKeyword('LABELS:', cs ) )    
             aux=cs(length('labels:  '):end);
             ind=find(aux~='"' & aux~=' '); 
             headerInfo.labels=aux(ind);
        elseif ( foundKeyword('SPACE:', cs ) )
            % Starts defining orientation (either space or space dimension, not
            % both)
            if isfield(headerInfo, 'spacedimension')
                fprintf('WARNING nhdr_nrrd_read %s: ignored field ''space'' because ''space dimension'' was specified before.\n', fopen(fidr));
            else
                headerInfo.space = strtrim( cs( length('SPACE:')+1:end ) );
                internal_spacedimension = nrrd_getSpaceDimensions(headerInfo.space);  % for internal use only
                assert(internal_spacedimension > 0, sprintf('nhdr_nrrd_read %s: unrecognized ''space'' descriptor ''%s''\n',fopen(fidr), headerInfo.space));
                define_orientation = 1;
            end
            
        elseif ( foundKeyword('SPACE DIMENSION:', cs) )
            % Starts defining orientation (either space or space dimension, not
            % both)
            if isfield(headerInfo, 'space')
                fprintf('WARNING nhdr_nrrd_read %s: ignored field ''space dimension'' because ''space'' was specified before.\n', fopen(fidr));
            else
                headerInfo.spacedimension = sscanf( cs( length('SPACE DIMENSION:')+1:end), '%i' );
                internal_spacedimension = headerInfo.spacedimension;  % for internal use only
                assert(numel(headerInfo.spacedimension)==1,...
                        sprintf('nhdr_nrrd_read %s: ''space dimension'' should be specified as one integer number. Found %d element(s) instead.\n',...
                                fopen(fidr), numel(headerInfo.spacedimension)));
                assert(headerInfo.spacedimension > 0, ...
                        sprintf('nhdr_nrrd_read %s: ''space dimension'' should be specified as a strictly positive integer (found %d).\n',...
                                fopen(fidr),headerInfo.spacedimension));
                define_orientation = 1;
            end
            
        elseif ( foundKeyword('SPACE DIRECTIONS:', cs ) )
            % Required if orientation defined but must come after space or space
            % dimension
            % space directions: <vector[0]> <vector[1]> ... <vector[dim-1]>
            
            assert(define_orientation==1, ...
                    sprintf('nhdr_nrrd_read %s: field ''space directions'' cannot be specified before ''space'' or ''space dimension''.\n',fopen(fidr)));
            % FIXME: ideally, check that the sum of elements including none and
            % (.,.,.) matches headerInfo.dimension, i.e. the total dimension.This
            % script currently gladly accepts absurd strings such as '1 0, 0 0 1 ( 0 none  0 ) 0 1,' 
            % instead of '(1,0,0) none none (0,1,0)  (0,0,1) none ' for instance. 
            % Maybe a problem ? Many people seem to mis-specify this
            % field so it may be a good idea to cut them some slack.
            space_dir_tmp = strtrim(cs(length('SPACE DIRECTIONS:')+1:end));     % remove leading and trailing spaces
            headerInfo.spacedirections = strsplit(space_dir_tmp);               % cell array of strings, ideally of the form {'(1,0,0)' '(0,2,0.1)' 'none'  '(0.1,0.1,1.1)'} if internal_spacedimension==3
            SD_data = strrep(space_dir_tmp, 'none', '' );                       % ignore "none" entries
            SD_data = extractNumbersWithout(SD_data, {'(',')',','} );           % detects numbers up to first non numerical entry
            assert(numel(SD_data) == (internal_spacedimension)^2,...
                   sprintf('nhdr_nrrd_read %s: expected ''space directions'' to specify a %d-by-%d matrix (%d total elements) in agreement with world space dimension. Found %d element(s) instead.\n',...
                            fopen(fidr),internal_spacedimension, internal_spacedimension, (internal_spacedimension)^2, numel(SD_data)) );
            headerInfo.spacedirections_matrix = reshape(SD_data(:), [internal_spacedimension, internal_spacedimension]);
            % correctness of field specification is checked below after 
            % whole header is parsed because we need to be sure that the 
            % dimension field was set
            
        elseif ( foundKeyword('SPACE UNITS:', cs ) )
            % Always optional, must come after space or space dimension
            
            assert(define_orientation==1,...
                    sprintf('nhdr_nrrd_read %s: field ''space units'' cannot be specified before ''space'' or ''space dimension''.\n', fopen(fidr)));
            
            space_units_tmp = strrep( cs(length('SPACE UNITS:')+1:end), 'none', '');        % ignore none entries
            % FIXME:  ideally, check that the sum of elements including none and
            % " " matches headerInfo.dimension, i.e. the total dimension, as
            % specified in the standard. Standard a bit unclear: should
            % unknown units be specified with "???", "none" or "" ? (quotes
            % seem to be required). 
            
            headerInfo.spaceunits = extract_spaceunits_list( space_units_tmp );
            
            assert(length(headerInfo.spaceunits) == internal_spacedimension,...
                    sprintf('nhdr_nrrd_read %s: expected ''space units'' to contain %d elements enclosed in double quotes to match the ''space'' or ''space dimension'' field but found the following %d element(s):\n%s\n',...
                            fopen(fidr), internal_spacedimension, length(headerInfo.spaceunits), sprintf('%s\t',headerInfo.spaceunits{:}) ));
            
        elseif ( foundKeyword('SPACE ORIGIN:', cs ) )
            % Always optional, must come after space or space dimension
            
            assert(define_orientation==1, sprintf('nhdr_nrrd_read %s: field ''space origin'' cannot be specified before ''space'' or ''space dimension''.\n', fopen(fidr)));
            
            iSO = extractNumbersWithout( cs(length('SPACE ORIGIN:')+1:end), {'(',')',','} );
            
            assert(numel(iSO) == internal_spacedimension,...
                    sprintf('nhdr_nrrd_read %s: expected ''space origin'' to specify a %d-element vector to match the ''space'' or ''space dimension'' field but found %d element(s).\n',...
                            fopen(fidr),internal_spacedimension, numel(iSO)) );
            
            headerInfo.spaceorigin = iSO(:);
            
        elseif ( foundKeyword('MEASUREMENT FRAME:', cs ) )
            % Always optional, must come after space or space dimension
            
            assert(define_orientation==1,...
                    sprintf('nhdr_nrrd_read %s: field ''measurement frame'' cannot be specified before ''space'' or ''space dimension''.\n', fopen(fidr)));
            
            measframe_str = strrep( cs(length('MEASUREMENT FRAME:')+1:end), 'none', '');
            
            iMF = extractNumbersWithout( measframe_str, {'(',')',','} ); % fails if non number entries are not previously removed
            
            assert(numel(iMF) == (internal_spacedimension)^2,...
                sprintf('nhdr_nrrd_read %s: expected ''measurement frame'' to specify a %d-by-%d matrix (%d total elements) but found %d element(s).\n',...
                        fopen(fidr), internal_spacedimension, internal_spacedimension, (internal_spacedimension)^2, numel(iMF)));
            
            headerInfo.measurementframe = reshape(iMF(:), [internal_spacedimension, internal_spacedimension]);
            
            % FIXME: this will gladly accept '1 0 0 0 1 0 0 0 1' instead of
            % '(1,0,0) (0,1,0)  (0,0,1)' for instance. But it should have
            % length spacedimension according to standard so this is not too bad.
            
        elseif ( foundKeyword('THICKNESSES:', cs ) )
            
            sThicknesses = extractStringList( cs(length('THICKNESSES:')+1:end) ); % fixed bug with extractStringList where 2 entries are present
            iThicknesses = [];
            lenThicknesses = length( sThicknesses );
            for iI=1:lenThicknesses
                iThicknesses = [iThicknesses, str2double( sThicknesses{iI} ) ];
            end
            headerInfo.thicknesses = iThicknesses;
            
        elseif ( foundKeyword('CENTERINGS:', cs ) )
            
            headerInfo.centerings = extractStringList( cs(length('CENTERINGS:')+1:end ) ); % fixed bug with extractStringList where 2 entries are present
            
        elseif ( foundKeyword('LINE SKIP:', cs) || foundKeyword('LINESKIP', cs) )
            
            if foundKeyword('LINE SKIP:', cs)
                headerInfo.lineskip = sscanf( cs( length('LINE SKIP:')+1:end ), '%d' );
            else
                headerInfo.lineskip = sscanf( cs( length('LINESKIP:')+1:end ), '%d' );
            end
            assert(headerInfo.lineskip >= 0,...
                    sprintf('Field lineskip or line skip should be greater than or equal to zero, detected %d.\n', headerInfo.lineskip));
        elseif ( foundKeyword('BYTE SKIP:', cs) || foundKeyword('BYTESKIP:', cs) )
            
            if foundKeyword('BYTE SKIP:', cs)
                headerInfo.byteskip = sscanf( cs( length('BYTE SKIP:')+1:end ), '%d' );
            else
                headerInfo.byteskip = sscanf( cs( length('BYTESKIP:')+1:end ), '%d' );
            end
            assert(headerInfo.byteskip >= -1, sprintf('Field byteskip or byte skip can only take non-negative integer values or -1, detected %d.\n',headerInfo.byteskip));
            
        elseif ( foundKeyword('MODALITY', cs ) )
            
            headerInfo.modality = strtrim( extractKeyValueString( cs(length('MODALITY')+1:end ) ) );
            
        elseif ( foundKeyword('DWMRI_B-VALUE', cs ) )
            
            headerInfo.bvalue = str2double( extractKeyValueString( cs(length('DWMRI_B-VALUE')+1:end ) ) );
            
        elseif ( foundKeyword('DWMRI_GRADIENT_', cs ) )
            
            [iGNr, dwiGradient] = extractGradient( cs(length('DWMRI_GRADIENT_')+1:end ) );
            headerInfo.gradients(iGNr+1,:) = dwiGradient;
            % FIXME: make it more general than numbering limited to 4 digits?
            
        elseif ( foundKeyword('DATA FILE:', cs ) || foundKeyword('DATAFILE:', cs) )
            % This tells us it is a detached header and ends it. 3 possibilities here
            
            field_value = strtrim( cs(length('DATA FILE:')+1:end) );
            if foundKeyword('DATAFILE:', cs)
                field_value = strtrim( cs(length('DATAFILE:')+1:end) );
            end
            
            [filelist, LIST_mode, subdim] = extract_datafiles(field_value);
            
            
            headerInfo.datafiles = filelist;
            
            % In LIST mode, filelist is empty and is filled by reading the
            % rest of the header
            if LIST_mode
                datafile_cnt = 0;
                while ~feof(fidr)
                    cs = fgetl(fidr);
                    if isempty(cs)
                        break;
                    end
                    datafile_cnt = datafile_cnt + 1;
                    headerInfo.datafiles{datafile_cnt} = cs;
                end
            end
            
            % We could technically break the loop here bc data file/datafile is
            % supposed to close a header; however I think it does no harm to
            % keep parsing (in LIST_mode, end of file reached anyways)
            
        else
            
            % see if we are dealing with a comment
            
            csTmp = strtrim( cs );
            if ( csTmp(1)~='#' && ~strcmp(cs(1:4),'NRRD') )
                fprintf('WARNING nhdr_nrrd_read: Could not parse input line: ''%s'' \n', cs );
                % REVIEW: is it better to blindly write it to the output
                % structure in order to be able to write the exact same file
                % later on ?
            end
        end
        
    end
    
    
    % Check for required fields 
    % REVIEW: should this be checked only when the data is read? People 
    % might be interested in just parsing the header no matter what is 
    % in it...)
    assert(isfield(headerInfo, 'sizes'), 'Missing required ''sizes'' field in header');
    assert(isfield(headerInfo, 'dimension'), 'Missing required ''dimension'' field in header');
    assert(isfield(headerInfo, 'encoding'), 'Missing required ''encoding'' field in header');
    assert(isfield(headerInfo, 'type'), 'Missing required ''type'' field in header');
    % Check other cross-field dependencies, etc. 
    % TODO: assert lengths of all
    % fields specified, check that kinds and axis sizes match, etc.
    if isfield(headerInfo, 'byteskip') && headerInfo.byteskip == -1
        assert( strcmpi(headerInfo.encoding, 'raw'), sprintf('byte skip value of -1 is only valid with raw encoding. See definition of NRRD File Format for more information.\n'));
    end
           
    matlabdatatype = nrrd_getMatlabDataType(headerInfo.type);
    
    %  endian field required if data defined on 2 bytes or more
    if ~(any(strcmpi(headerInfo.encoding,{'txt', 'text', 'ascii'})) || any(strcmpi(matlabdatatype, {'int8', 'uint8'})))
        assert(isfield(headerInfo, 'endian'), 'Missing required ''endian'' field in header');
    else
        if ~isfield(headerInfo,'endian')
            headerInfo.endian = 'little';    % useless but harmless, just for code compatibility because endian field may be accessed later on while reading data
        end
    end
    
    % Space and orientation information
    if define_orientation
        assert(isfield(headerInfo, 'spacedirections'), 'Missing field ''space directions'', required if either ''space'' or ''space dimension'' is set.');
        if length(headerInfo.spacedirections) ~= headerInfo.dimension
            fprintf('WARNING nhdr_nrrd_read %s:\n Unexpected format found for ''space directions'' specifier.\n Expected a combination of none and vectors delimited by parentheses such as (1.2,0.2,0.25), see definition of NRRD File Format for more information.\n', fopen(fidr));
        end
    end
    
    % TODO: add support for positive line skips 
    if isfield(headerInfo, 'lineskip') && headerInfo.lineskip > 0
        assert(isfield(headerInfo, 'byteskip') && headerInfo.byteskip == -1, sprintf('lineskip option is currently not supported and can only be set to zero unless raw encoding is used and byte skip is set to -1 (which cancels the effect of lineskip altogether).\n'));
    end  
    
    % TODO: add support for positive byte skips (see line skip)
    if isfield(headerInfo, 'byteskip') && ~strcmpi(headerInfo.encoding,'raw')
        assert( headerInfo.byteskip == 0, sprintf('byte skip option with non raw encoding is currently not supported and can only be set to zero.\n'));
    end
    if isfield(headerInfo, 'byteskip') && strcmpi(headerInfo.encoding, 'raw')
        assert( headerInfo.byteskip == -1, sprintf('non-negative byte skip values with raw encoding are currently not supported; byte skip can only be set to -1.\n'));
    end
    
catch me
    % Clean up before raising error
    fclose(fidr);
    rethrow(me);
end
% Read the data. Detect if data is in the file or in a detached data file
% and check what type of detached file it is if need be.
if bReadData
    N_data_tot = prod(headerInfo.sizes);
    if isfield(headerInfo, 'datafiles')
        % This is a detached header file
        
        % We no longer need to read the header (closing it before reading
        % all the detached data files is safer)
        fclose(fidr);           
        
        % TODO: we currently only read slices along the slowest axis (i.e.,
        % last coordinate)
        if ~isempty(subdim) && subdim~=headerInfo.dimension
            error('nhdr_nrrd_read (detached header): reading data from slices along axis other than the last (i.e. slowest) one, is currently not supported.\nLast argument [<subdim>] in ''data file'' or ''datafile'' field should be removed or set equal to %d, which is the detected ''dimension'' field value, for now.\n',...
                        headerInfo.dimension);
        end
        
        % Read data chunk by chunk from detached data files
        N_data_files = length(headerInfo.datafiles);
        assert(mod(N_data_tot, N_data_files)==0, sprintf('Number of detected data files (%d) does not divide total number of values contained in data %d obtained from prod(sizes=[%s]).\n',...
                                                                        N_data_files, N_data_tot, sprintf('%d ',headerInfo.sizes)));
        N_data_per_file = N_data_tot/N_data_files;
        
        headerInfo.data = zeros(headerInfo.sizes, matlabdatatype);      % specify right type of zeros, otherwise double by default
        
        for i = 1:N_data_files
            
            % Check type of detached data file
            [~,fname_data,ext_data] = fileparts(headerInfo.datafiles{i});        % redundant because done above
            
            data_ind = (i-1)*N_data_per_file+1:i*N_data_per_file;
            
            if strcmpi(ext_data,'.nhdr')
                
                error('datafile %di/%d: nhdr file should not be used as detached data file.\n',i,length(headerInfo.datafiles));
                
            elseif strcmpi(ext_data, '.nrrd')
                % Detached NRRD file
                
                bRead_Detached_Data = true;
                tmp_struct = nhdr_nrrd_read(fullfile(mainFpath, [fname_data, ext_data]), bRead_Detached_Data);   % recursive call
                % TODO: check the rest of the structure for
                % inconsistencies with metadata from header
                assert(N_data_files==1 || headerInfo.dimension == tmp_struct.dimension +1, sprintf(...
                    'nhdr_nrrd_read detached header %s: the number of dimensions in detached nrrd data file %d/%d (%s) should be one fewer than that of the header file.\nDetected %d instead of %d.\nDifferent datafile dimensions as specified by the nrrd standard are not supported as of now. \n',...
                    fopen(fidr), i, N_data_files, headerInfo.datafiles{i}, tmp_struct.dimension, headerInfo.dimension-1 ));
                
                headerInfo.data(data_ind) = tmp_struct.data(:);
                
                % Store detached data header for last detached file seen,
                % assuming that all are the same. Do not store data though
                % as it would be redundant.
                tmp_struct = rmfield(tmp_struct, 'data');
                headerInfo.detached_header = tmp_struct;
            else
                % e.g., detached .raw file
                
                fid_data = fopen( fullfile(mainFpath, [fname_data, ext_data]), 'r');
                if( fid_data < 1 )
                    error('While reading detached header file %s: detached data file number %d/%d (%s) could not be opened.\n', nrrdFileName, i, N_data_files, headerInfo.datafiles{i});
                end
                try
                    tmp_data = readData(fid_data, N_data_per_file, headerInfo.encoding, matlabdatatype);
                    fclose(fid_data);
                catch me_detached
                    fclose(fid_data);
                    rethrow(me_detached);
                end
                
                tmp_data = adjustEndian(tmp_data, headerInfo.endian);
                
                headerInfo.data(data_ind) = tmp_data(:);
            end
            
        end
        
    else
        % This is a NRRD standalone file: read data directly from it
        try
            headerInfo.data = readData(fidr, N_data_tot, headerInfo.encoding, matlabdatatype);
            fclose(fidr);
        catch me_detached
            fclose(fidr);
            rethrow(me_detached);
        end        
        headerInfo.data = adjustEndian(headerInfo.data, headerInfo.endian);
        headerInfo.data = reshape(headerInfo.data, headerInfo.sizes(:)');   % data into expected form. Transpose required by reshape function bc size vectors need to be row vectors
    end
else
    fclose(fidr);
end
end
% -------------------------------------------------------------------%
% ====================================================================
% --- Auxiliary functions -------------------------------------------%
% ====================================================================
function [iGNr, dwiGradient] = extractGradient( st )
% first get the gradient number
iGNr = str2num( st(1:4) );           % FIX numbering limited to 4 digits (synchronize with writer module)
% find where the assignment is
assgnLoc = strfind( st, ':=' );
if ( isempty(assgnLoc) )
    dwiGradient = [];
    return;
else
    
    dwiGradient = sscanf( st(assgnLoc+2:end), '%f' );
    
end
end
% Return part of string after :=
function kvs = extractKeyValueString( st )
assgnLoc = strfind( st, ':=' );
if ( isempty(assgnLoc) )
    kvs = [];
    return;
else
    
    kvs = st(assgnLoc(1)+2:end);
    
end
end
% Turn space-separated list into cell array of strings. 
function sl = extractStringList( strList )
sl = strsplit(strtrim(strList)); % old Matlab file exchange version had a strange bug with lists of length 2
end
% Store in an array the list of numbers separated by the tokens listed in
% the withoutTokens cell array 
function iNrs = extractNumbersWithout( inputString, withoutTokens )
auxStr = inputString;
for iI=1:length( withoutTokens )
    
    auxStr = strrep( auxStr, withoutTokens{iI}, ' ' );
    
end
iNrs = sscanf( auxStr, '%f' );
end
% Return true if the string keyWord is the beginning of the content string
% cs (ignoring case)
function fk = foundKeyword( keyWord, cs )
lenKeyword = length( keyWord );
fk = (lenKeyword <= length(cs)) && strcmpi( cs(1:lenKeyword), keyWord);
end
% Return cell array of strings with space units in the form {mm, km, m}
% from input of the form "mm " "mm" " m" where undesired blank spaces may
% have crept in. The double quotes are removed in the process but will be
% added by the nhdr/nrrd writer function.
function su_ca = extract_spaceunits_list( fieldValue )
fv_trimmed = strtrim( fieldValue );
su_ca = strsplit(fv_trimmed, '"');                              % units are delimited by double quotes
su_ca = su_ca(~ ( strcmp(su_ca, '') | strcmp(su_ca, ' ') ) );   % remove empty or blank space strings
for i = 1:length(su_ca)
    su_ca{i} = strtrim( su_ca{i} );
end
end
% Extract data file list into a cell array from the datafile or data file
% field in a NHDR file, stripped off its leading and trailing spaces;
% LIST_mode is set to one if a list of files closes the header;
% subdim is empty by default and may be set through either of these field
% specifications:
% data file: <format> <min> <max> <step> [<subdim>]
% data file: LIST [<subdim>]
function [filelist, LIST_mode, subdim] = extract_datafiles( field_string)
field_string = strtrim(field_string);
filelist = {};
LIST_mode = 0;
subdim = [];
if length(field_string)>= 4 && strcmpi(field_string(1:4), 'LIST')
    % LIST mode
    subdim = sscanf(field_string(5:end), '%d'); % assert 1<=dimensions 
    LIST_mode = 1;
    return;
else
    % single detached data file or multiple files written in concise form, non LIST mode
    str_lst = strsplit(field_string); % without delimiters specified, splits at any sequence in the set  {' ','\f','\n','\r','\t','\v'}
    if length(str_lst) == 1
        % Single detached data file:
        filelist{1} = str_lst{1};
        % subdim does not make sense here since all of the data is
        % contained in the deatached data file
    elseif any(length(str_lst)==[0, 2, 3]) || length(str_lst)>5 
        % Invalid cases (2, 3 or striclty more than 5 entries)
        error('error in ''data list'' (or ''data list'') field, in non ''LIST'' mode:\nexpected  ''<filename>'' or ''<format> <min> <max> <step> [<subdim>]'' but found instead ''%s'' (%d elements instead of 1, 4 or 5).\n',...
            sprintf('%s ',str_lst{:}), length(str_lst));
    else 
        % Multiple detached data files with filenames including integral
        % values generated by min:step:max
                
        str_format = str_lst{1};
        id_min = sscanf(str_lst{2}, '%d');
        id_max = sscanf(str_lst{3}, '%d');
        step = sscanf(str_lst{4}, '%d');
        
        % check format and indexing values
        assert(step~=0,sprintf('detached data files with names specified by sprintf()-like format:  step should be strictly positive or negative, not zero.\n'));
        if step < 0
            assert(id_min >= id_max, sprintf('detached data files with names specified by sprintf()-like format:  when step is <0, min should be larger than or equal to max, here we found min=%d < max=%d.\n',id_min, id_max));
        else
            assert(id_min <= id_max, sprintf('detached data files with names specified by sprintf()-like format: : when step is >0, min should be smaller than or equal to max, here we found min=%d > max=%d.\n',id_min, id_max));
        end
         
        % populate data file list
        fileIDs = id_min:step:id_max;
        filelist = cell(length(fileIDs), 1);
        for i = 1:length(fileIDs)  % does not necessarily incude <max>, it cannot be larger than <max> 
            expanded_fname = sprintf(str_format, fileIDs(i));
            filelist{i} = expanded_fname;
        end
        
        if length(str_lst) == 5
            subdim = sscanf(str_lst{5}, '%d');
        end   
    end
end
        
end
% Read data in a column vector from open file with pointer fidIn at the
% beginning of the data to read. The other arguments specify the number of
% elements to read (determined in main body of reader function above), data
% file encoding (specified in NHDR/NRRD header) and the Matlab type the
% data should be converted to (determined in main body of reader function
% above).
function data = readData(fidIn, Nelems, meta_encoding, matlabdatatype)
switch (meta_encoding)
    case {'raw'}
        % This implicitely assumes that byteskip was set to -1 in raw
        % encoding
        data_bytes = (Nelems)*sizeOf(matlabdatatype);         % bytes of useful info to return in data array
        fseek(fidIn,0,'eof');                           % set position pointer to end of file
        tot_bytes_file = ftell(fidIn);                  % current location of the position pointer (I believe as a byte offset from beginning of file)
        fseek(fidIn,tot_bytes_file-data_bytes,'bof');   % make sure position pointer is right before useful data starts
        
        % Just read binary file (original version of the code)
        data = fread(fidIn, inf, [matlabdatatype '=>' matlabdatatype]);
        
    case {'gzip', 'gz'}
        
        tmp = fread(fidIn, Inf, 'uint8=>uint8'); % byte per byte
        
        tmpBase = tempname(pwd);
        tmpFile = [tmpBase '.gz'];
        fidTmp = fopen(tmpFile, 'wb');  % this creates tmpFile
        assert(fidTmp > 3, 'Could not open temporary file for GZIP decompression')
        
        try
            fwrite(fidTmp, tmp, 'uint8');
        catch me
            fclose(fidTmp);
            delete(tmpFile);
            rethrow(me);
        end
        fclose(fidTmp);        
        
        try
            gunzip(tmpFile); % this creates tmpBase (tmpFile without the .gz)
        catch me
            delete(tmpFile);
            rethrow(me);
        end
        delete(tmpFile);
        
        fidTmp = fopen(tmpBase, 'rb');
        % cleaner = onCleanup(@() fclose(fidTmp));
        try
            data = readData(fidTmp, Nelems, 'raw', matlabdatatype);    % recursive call            
        catch me
            fclose(fidTmp);
            delete(tmpBase);
            rethrow(me);
        end
        fclose(fidTmp);
        delete(tmpBase);
    case {'txt', 'text', 'ascii'}
        
        data = fscanf(fidIn, '%f');
        data = cast(data, matlabdatatype);
        
    otherwise
        assert(false, 'Unsupported encoding')
end
% Check number of read elements (sometimes there is an offest of about 1)
assert(Nelems == numel(data), sprintf('Error reading binary content of %s: detected %d elements instead of %d announced in header.\n',fopen(fidIn), numel(data), Nelems));
end
% Switch byte ordering according to endianness specified in nhdr/nrrd file.
% The metadata should just contain a field 'endian' set to 'little' or
% 'big'.
function data = adjustEndian(data, meta_endian)
[~,~,endian] = computer();
needToSwap = (isequal(endian, 'B') && isequal(lower(meta_endian), 'little')) || ...
    (isequal(endian, 'L') && isequal(lower(meta_endian), 'big'));
if (needToSwap)
    data = swapbytes(data);
end
end
% Size of Matlab's numeric classes in bytes
% See: https://stackoverflow.com/questions/16104423/size-of-a-numeric-class
function numBytes = sizeOf(dataClass)
    % create temporary variable of data type specified
    eval(['var = ' dataClass '(0);']); 
    % Use the functional form of whos, and get the number of bytes   
    W = whos('var');
    numBytes = W.bytes;
end
% Get matlab type corresponding to each nrrd type
% Date: November 2017
% Author: Gaetan Rensonnet
function matlabdatatype = nrrd_getMatlabDataType(metaType)
% Determine the datatype
switch (metaType)
    case {'signed char', 'int8', 'int8_t'}
        matlabdatatype = 'int8';
        
    case {'uchar', 'unsigned char', 'uint8', 'uint8_t'}
        matlabdatatype = 'uint8';
        
    case {'short', 'short int', 'signed short', 'signed short int', ...
            'int16', 'int16_t'}
        matlabdatatype = 'int16';
        
    case {'ushort', 'unsigned short', 'unsigned short int', 'uint16', ...
            'uint16_t'}
        matlabdatatype = 'uint16';
        
    case {'int', 'signed int', 'int32', 'int32_t'}
        matlabdatatype = 'int32';
        
    case {'uint', 'unsigned int', 'uint32', 'uint32_t'}
        matlabdatatype = 'uint32';
        
    case {'longlong', 'long long', 'long long int', 'signed long long', ...
            'signed long long int', 'int64', 'int64_t'}
        matlabdatatype = 'int64';
        
    case {'ulonglong', 'unsigned long long', 'unsigned long long int', ...
            'uint64', 'uint64_t'}
        matlabdatatype = 'uint64';
        
    case {'float'}
        matlabdatatype = 'single';
        
    case {'double'}
        matlabdatatype = 'double';
        
    case {'block'}
        error('Data type ''block'' is currently not supported.\n');
        
    otherwise
        assert(false, 'Unknown datatype')
end
end
% Get number of space dimensions ("space dimension" field) according to the
% space descriptor, following the format definition at
% http://teem.sourceforge.net/nrrd/format.html#spacedirections.
%
% Returns 3 if "space" field is one of
% "right-anterior-superior" or "RAS"
% "left-anterior-superior" or "LAS"
% "left-posterior-superior" or "LPS"
% "scanner-xyz"
% "3D-right-handed"
% "3D-left-handed"
%
% Returns 4 if "space" field is one of
% "right-anterior-superior-time" or "RAST"
% "left-anterior-superior-time" or "LAST"
% "left-posterior-superior-time" or "LPST"
% "scanner-xyz-time"
% "3D-right-handed-time"
% "3D-left-handed-time"
%
% Date: October 25, 2017
% Author: Gaetan Rensonnet
function sd = nrrd_getSpaceDimensions(spacedescriptor)
if any(strcmpi(spacedescriptor,...
        {'right-anterior-superior', 'RAS',...
        'left-anterior-superior', 'LAS',...
        'left-posterior-superior', 'LPS',...
        'scanner-xyz',...
        '3D-right-handed',...
        '3D-left-handed'}...
        ))
    sd = 3;
    
elseif any(strcmpi(spacedescriptor,...
        {'right-anterior-superior-time', 'RAST',...
        'left-anterior-superior-time', 'LAST',...
        'left-posterior-superior-time', 'LPST', ...
        'scanner-xyz-time', ...
        '3D-right-handed-time', ...
        '3D-left-handed-time'}...
        ))
    sd = 4;
    
else
    sd = -1;
    % Unrecognized nrrd space descriptor (grace under fire)
end
end
