function nii = mcxloadnii(filename)
%
%    nii=mcxloadnii(filename)
%
%    author: Qianqian Fang (q.fang <at> neu.edu)
%
%    input:
%        fname: the file name to the .nii file
%
%    output:
%        nii.img: the data volume read from the nii file
%        nii.datatype: the data type of the voxel, in matlab data type string
%        nii.datalen: data count per voxel - for example RGB data has 3x
%                    uint8 per voxel, so datatype='uint8', datalen=3
%        nii.voxelbyte: total number of bytes per voxel: for RGB data,
%                    voxelbyte=3; also voxelbyte=header.bitpix/8
%        nii.hdr: file header info, a structure has the full nii header
%                    key subfileds include
%
%              sizeof_hdr: must be 348
%              dim: short array, dim(2: dim(1)+1) defines the array size
%              datatype: the type of data stored in each voxel
%              bitpix: total bits per voxel
%              magic: must be 'ni1\0' or 'n+1\0'
%
%              For the detailed nii header, please see
%              https://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1.h
%
%    this file is part of Monte Carlo eXtreme (MCX)
%    License: GPLv3, see https://mcx.space for details
%

header = memmapfile(filename,                ...
                    'Offset', 0,                           ...
                    'Writable', false,                     ...
                    'Format', {                            ...
                               'int32'   [1 1]  'sizeof_hdr'      % !< MUST be 348           %  % int sizeof_hdr;       %  ...
                               'int8'    [1 10] 'data_type'       % !< ++UNUSED++            %  % char data_type[10];   %  ...
                               'int8'    [1 18] 'db_name'         % !< ++UNUSED++            %  % char db_name[18];     %  ...
                               'int32'   [1 1]  'extents'         % !< ++UNUSED++            %  % int extents;          %  ...
                               'int16'   [1 1]  'session_error'   % !< ++UNUSED++            %  % short session_error;  %  ...
                               'int8'    [1 1]  'regular'         % !< ++UNUSED++            %  % char regular;         %  ...
                               'int8'    [1 1]  'dim_info'        % !< MRI slice ordering.   %  % char hkey_un0;        %  ...
                               'uint16'  [1 8]  'dim'             % !< Data array dimensions.%  % short dim[8];         %  ...
                               'single'  [1 1]  'intent_p1'       % !< 1st intent parameter. %  % short unused8/9;      %  ...
                               'single'  [1 1]  'intent_p2'       % !< 2nd intent parameter. %  % short unused10/11;    %  ...
                               'single'  [1 1]  'intent_p3'       % !< 3rd intent parameter. %  % short unused12/13;    %  ...
                               'int16'   [1 1]  'intent_code'     % !< NIFTI_INTENT_* code.  %  % short unused14;       %  ...
                               'int16'   [1 1]  'datatype'        % !< Defines data type!    %  % short datatype;       %  ...
                               'int16'   [1 1]  'bitpix'          % !< Number bits/voxel.    %  % short bitpix;         %  ...
                               'int16'   [1 1]  'slice_start'     % !< First slice index.    %  % short dim_un0;        %  ...
                               'single'  [1 8]  'pixdim'          % !< Grid spacings.        %  % float pixdim[8];      %  ...
                               'single'  [1 1]  'vox_offset'      % !< Offset into .nii file %  % float vox_offset;     %  ...
                               'single'  [1 1]  'scl_slope'       % !< Data scaling: slope.  %  % float funused1;       %  ...
                               'single'  [1 1]  'scl_inter'       % !< Data scaling: offset. %  % float funused2;       %  ...
                               'int16'   [1 1]  'slice_end'       % !< Last slice index.     %  % float funused3;       %  ...
                               'int8'    [1 1]  'slice_code'      % !< Slice timing order.   %                             ...
                               'int8'    [1 1]  'xyzt_units'      % !< Units of pixdim[1..4] %                             ...
                               'single'  [1 1]  'cal_max'         % !< Max display intensity %  % float cal_max;        %  ...
                               'single'  [1 1]  'cal_min'         % !< Min display intensity %  % float cal_min;        %  ...
                               'single'  [1 1]  'slice_duration'  % !< Time for 1 slice.     %  % float compressed;     %  ...
                               'single'  [1 1]  'toffset'         % !< Time axis shift.      %  % float verified;       %  ...
                               'int32'   [1 1]  'glmax'           % !< ++UNUSED++            %  % int glmax;            %  ...
                               'int32'   [1 1]  'glmin'           % !< ++UNUSED++            %  % int glmin;            %  ...
                               'int8'    [1 80] 'descrip'         % !< any text you like.    %  % char descrip[80];     %  ...
                               'int8'    [1 24] 'aux_file'        % !< auxiliary filename.   %  % char aux_file[24];    %  ...
                               'int16'   [1 1]  'qform_code'      % !< NIFTI_XFORM_* code.   %  %-- all ANALYZE 7.5 --- %  ...
                               'int16'   [1 1]  'sform_code'      % !< NIFTI_XFORM_* code.   %  %below here are replaced%  ...
                               'single'  [1 1]  'quatern_b'       % !< Quaternion b param.   %             ...
                               'single'  [1 1]  'quatern_c'       % !< Quaternion c param.   %             ...
                               'single'  [1 1]  'quatern_d'       % !< Quaternion d param.   %             ...
                               'single'  [1 1]  'qoffset_x'       % !< Quaternion x shift.   %             ...
                               'single'  [1 1]  'qoffset_y'       % !< Quaternion y shift.   %             ...
                               'single'  [1 1]  'qoffset_z'       % !< Quaternion z shift.   %             ...
                               'single'  [1 4]  'srow_x'          % !< 1st row affine transform.   %           ...
                               'single'  [1 4]  'srow_y'          % !< 2nd row affine transform.   %           ...
                               'single'  [1 4]  'srow_z'          % !< 3rd row affine transform.   %           ...
                               'int8'    [1 16] 'intent_name'     % !< 'name' or meaning of data.  %           ...
                               'int8'    [1 4]  'magic'           % !< MUST be "ni1\0" or "n+1\0". %           ...
                               'int8'    [1 4]  'extension'       % !< header extension      %             ...
                              });

nii.hdr = header.Data(1);

type2byte = [
             0  0  % unknown                      %
             1  0  % binary (1 bit/voxel)         %
             2  1  % unsigned char (8 bits/voxel) %
             4  2  % signed short (16 bits/voxel) %
             8  4  % signed int (32 bits/voxel)   %
             16  4  % float (32 bits/voxel)        %
             32  8  % complex (64 bits/voxel)      %
             64  8  % double (64 bits/voxel)       %
             128  3  % RGB triple (24 bits/voxel)   %
             255  0  % not very useful (?)          %
             256  1  % signed char (8 bits)         %
             512  2  % unsigned short (16 bits)     %
             768  4  % unsigned int (32 bits)       %
             1024  8  % long long (64 bits)          %
             1280  8  % unsigned long long (64 bits) %
             1536 16  % long double (128 bits)       %
             1792 16  % double pair (128 bits)       %
             2048 32  % long double pair (256 bits)  %
             2304  4  % 4 byte RGBA (32 bits/voxel)  %
            ];

type2str = {
            'uint8'    0   % unknown                       %
            'uint8'    0   % binary (1 bit/voxel)          %
            'uint8'    1   % unsigned char (8 bits/voxel)  %
            'uint16'   1   % signed short (16 bits/voxel)  %
            'int32'    1   % signed int (32 bits/voxel)    %
            'float32'  1   % float (32 bits/voxel)         %
            'float32'  2   % complex (64 bits/voxel)       %
            'float64'  1   % double (64 bits/voxel)        %
            'uint8'    3   % RGB triple (24 bits/voxel)    %
            'uint8'    0   % not very useful (?)           %
            'int8'     1   % signed char (8 bits)          %
            'uint16'   1   % unsigned short (16 bits)      %
            'uint32'   1   % unsigned int (32 bits)        %
            'long'     1   % long long (64 bits)           %
            'ulong'    1   % unsigned long long (64 bits)  %
            'uint8'    16  % long double (128 bits)        %
            'uint8'    16  % double pair (128 bits)        %
            'uint8'    32  % long double pair (256 bits)   %
            'uint8'    4   % 4 byte RGBA (32 bits/voxel)   %
           };

typeidx = find(type2byte(:, 1) == nii.hdr.datatype);

nii.datatype = type2str{typeidx, 1};
nii.datalen = type2str{typeidx, 2};
nii.voxelbyte = type2byte(typeidx, 2);

if (type2byte(typeidx, 2) == 0)
    nii.img = [];
    return
end

if (type2str{typeidx, 2} > 1)
    nii.hdr.dim = [nii.hdr.dim(1) + 1 uint16(nii.datalen) nii.hdr.dim(2:end)];
end

fid = fopen(filename, 'rb');
fseek(fid, nii.hdr.vox_offset, 'bof');
nii.img = fread(fid, prod(nii.hdr.dim(2:nii.hdr.dim(1) + 1)), [nii.datatype '=>' nii.datatype]);
fclose(fid);

nii.img = reshape(nii.img, nii.hdr.dim(2:nii.hdr.dim(1) + 1));
