From 2e4358648720de52f3f38446bc5fdcdfdbeef2d7 Mon Sep 17 00:00:00 2001 From: Qianqian Fang Date: Wed, 27 Mar 2024 23:02:03 -0400 Subject: [PATCH] [feat] merge nii/jnii/hdf5/tsv reading functions self-contained --- getvarfrom.m | 33 +++++ loadbidstsv.m | 92 ++++++++++++++ loadh5.m | 312 +++++++++++++++++++++++++++++++++++++++++++++++ loadjnifti.m | 57 +++++++++ loadnifti.m | 21 ++++ memmapstream.m | 81 ++++++++++++ nii2jnii.m | 306 ++++++++++++++++++++++++++++++++++++++++++++++ niicodemap.m | 93 ++++++++++++++ niiformat.m | 128 +++++++++++++++++++ niiheader2jnii.m | 65 ++++++++++ regrouph5.m | 134 ++++++++++++++++++++ transposemat.m | 40 ++++++ 12 files changed, 1362 insertions(+) create mode 100644 getvarfrom.m create mode 100644 loadbidstsv.m create mode 100644 loadh5.m create mode 100644 loadjnifti.m create mode 100644 loadnifti.m create mode 100644 memmapstream.m create mode 100644 nii2jnii.m create mode 100644 niicodemap.m create mode 100644 niiformat.m create mode 100644 niiheader2jnii.m create mode 100644 regrouph5.m create mode 100644 transposemat.m diff --git a/getvarfrom.m b/getvarfrom.m new file mode 100644 index 0000000..f0fbb60 --- /dev/null +++ b/getvarfrom.m @@ -0,0 +1,33 @@ +function p = getvarfrom(ws, name) +% +% p=getvarfrom(ws,name) +% +% get variable value by name from specified work-space +% +% author: Qianqian Fang, +% +% input: +% ws: name of the work-space, for example, 'base' +% name: name string of the variable +% +% output: +% p: the value of the specified variable, if the variable does not +% exist, return empty array +% +% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net) +% + +wsname = ws; +if (~iscell(ws)) + wsname = cell(1); + wsname{1} = ws; +end + +p = []; +for i = 1:length(wsname) + isdefined = evalin(wsname{i}, ['exist(''' name ''')']); + if (isdefined == 1) + p = evalin(wsname{i}, name); + break + end +end diff --git a/loadbidstsv.m b/loadbidstsv.m new file mode 100644 index 0000000..0eb6a10 --- /dev/null +++ b/loadbidstsv.m @@ -0,0 +1,92 @@ +function data = loadbidstsv(tsvfile, delim) +% +% data = loadbidstsv(tsvfile) +% or +% data = loadbidstsv(tsvfile, delim) +% +% Loading a BIDS-formatted .tsv (tab-separated values) or .tsv.gz file as a +% struct; numerical fields are converted to floating-point data records +% when possible; the header of the file is parsed to define sub-field +% names +% +% author: Qianqian Fang (q.fang neu.edu) +% +% input: +% tsvfile: the path to the .tsv file +% delim: (optional) if not set, tab ('\t') is used as column delimiter +% +% examples: +% data = loadbidstsv('participants.tsv'); +% +% license: +% BSD license, see LICENSE_BSD.txt files for details +% +% -- this function is part of JBIDS toolbox (https://neurojson.org/#software) +% + +if (nargin < 2) + delim = sprintf('\t'); +end + +data = struct; + +if (~isempty(regexp(tsvfile, '\.[Gg][Zz]$', 'once'))) + finput = fopen(tsvfile, 'rb'); + tsvdata = fread(finput, inf, 'uint8=>uint8'); + fclose(finput); + + if (~exist('gzipdecode', 'file')) + error('To process zipped files, you must install gzipdecode.m from the JSONLab toolbox: http://github.com/fangq/jsonlab'); + end + fid = char(gzipdecode(tsvdata)); + clear tsvdata; + [header, endpos] = regexp(fid, '([^\n\r]*)', 'once', 'tokens', 'end'); + if (~isempty(header)) + header = header{1}; + fid = fid((endpos + 1):end); + end +else + fid = fopen(tsvfile, 'rt'); + header = fgetl(fid); + header = regexprep(header, '\s*$', ''); +end + +if (isempty(header)) + return +end + +if (exist('strsplit')) + cols = strsplit(header, delim); +else + cols = regexp(header, '\t*([^\t]*)\t*', 'tokens'); + cols = cellfun(@(x) x{:}, cols, 'uniformoutput', 0); +end +cols = cellfun(@encodevarname, cols, 'uniformoutput', 0); +if (~isempty(cols)) + body = textscan(fid, [repmat('%s\t', 1, length(cols) - 1), '%s'], 'delimiter', '\t'); + if (length(body) ~= length(cols)) + error('invalid tsv'); + end + for i = 1:length(body) + try + % bodynum = cellfun(@(x) sscanf(regexprep(x,'^n/a$','NaN'), '%f'), body{i}, 'uniformoutput', 0); + bodynum = cellfun(@(x) sscanf(x, '%f\t'), body{i}, 'uniformoutput', 0); + if (exist('isna')) + len = cellfun(@(x) numel(x) * (~isna(sum(x))), bodynum); + else + len = cellfun(@numel, bodynum); + end + if (any(len)) + body{i}(len > 0) = bodynum(len > 0); + if (all(len)) + body{i} = cell2mat(body{i}); + end + end + catch ME + warning(ME.message); + end + data.(cols{i}) = body{i}(:).'; + end +end + +fclose(fid); diff --git a/loadh5.m b/loadh5.m new file mode 100644 index 0000000..048b07c --- /dev/null +++ b/loadh5.m @@ -0,0 +1,312 @@ +function varargout = loadh5(filename, varargin) +% +% [data, meta] = loadh5(filename) +% [data, meta] = loadh5(root_id) +% [data, meta] = loadh5(filename, rootpath) +% [data, meta] = loadh5(filename, rootpath, options) +% [data, meta] = loadh5(filename, options) +% [data, meta] = loadh5(filename, 'Param1',value1, 'Param2',value2,...) +% +% Load data in an HDF5 file to a MATLAB structure. +% +% author: Qianqian Fang (q.fang neu.edu) +% +% input +% filename +% Name of the file to load data from +% root_id: an HDF5 handle (of type 'H5ML.id' in MATLAB) +% rootpath : (optional) +% Root path to read part of the HDF5 file to load +% options: (optional) a struct or Param/value pairs for user specified options +% Order: 'creation' - creation order (default), or 'alphabet' - alphabetic +% Regroup: [0|1]: if 1, call regrouph5() to combine indexed +% groups into a cell array +% PackHex: [1|0]: convert invalid characters in the group/dataset +% names to 0x[hex code] by calling encodevarname.m; +% if set to 0, call getvarname +% ComplexFormat: {'realKey','imagKey'}: use 'realKey' and 'imagKey' +% as possible keywords for the real and the imaginary part +% of a complex array, respectively (sparse arrays not supported); +% a common list of keypairs is used even without this option +% Transpose: [1|0] - if set to 1 (default), the row-majored HDF5 +% datasets are transposed (to column-major) so that the +% output MATLAB array has the same dimensions as in the +% HDF5 dataset header. +% +% output +% data: a structure (array) or cell (array) +% meta: optional output to store the attributes stored in the file +% +% example: +% a={rand(2), struct('va',1,'vb','string'), 1+2i}; +% saveh5(a,'test.h5'); +% a2=loadh5('test.h5') +% a3=loadh5('test.h5','regroup',1) +% isequaln(a,a3.a) +% a4=loadh5('test.h5','/a1') +% +% This function was adapted from h5load.m by Pauli Virtanen +% This file is part of EasyH5 Toolbox: https://github.com/NeuroJSON/easyh5 +% +% License: GPLv3 or 3-clause BSD license, see https://github.com/NeuroJSON/easyh5 for details +% + +path = ''; +opt = struct; +if (bitand(length(varargin), 1) == 0) + opt = varargin2struct(varargin{:}); +elseif (length(varargin) >= 3) + path = varargin{1}; + opt = varargin2struct(varargin{2:end}); +elseif (length(varargin) == 1) + path = varargin{1}; +end + +opt.dotranspose = jsonopt('Transpose', 1, opt); +opt.stringarray = jsonopt('StringArray', 0, opt); +opt.rootpath = path; + +if (exist('OCTAVE_VERSION', 'builtin') ~= 0) + [varargout{1:nargout}] = load(filename, '-hdf5'); + if (opt.dotranspose) + varargout{1} = transposemat(varargout{1}); + end + return +end + +if (isa(filename, 'H5ML.id')) + loc = filename; +else + try + loc = H5F.open(filename, 'H5F_ACC_RDONLY', 'H5P_DEFAULT'); + catch ME + error('fail to open file'); + end +end + +if (~(isfield(opt, 'complexformat') && iscellstr(opt.complexformat) && numel(opt.complexformat) == 2)) + opt.complexformat = {'Real', 'Imag'}; +end + +opt.releaseid = 0; +vers = ver('MATLAB'); +if (~isempty(vers)) + opt.releaseid = datenum(vers(1).Date); +end + +if ((isfield(opt, 'order') && strcmpi(opt.order, 'alphabet')) || opt.releaseid < datenum('1-Jan-2015')) + opt.order = 'H5_INDEX_NAME'; +else + opt.order = 'H5_INDEX_CRT_ORDER'; +end + +try + if (nargin > 1 && ~isempty(path)) + try + rootgid = H5G.open(loc, path); + [varargout{1:nargout}] = load_one(rootgid, opt); + H5G.close(rootgid); + catch + [gname, dname] = fileparts(path); + rootgid = H5G.open(loc, gname); + [status, res] = group_iterate(rootgid, dname, struct('data', struct, 'meta', struct, 'opt', opt)); + if (nargout > 0) + varargout{1} = res.data; + elseif (nargout > 1) + varargout{2} = res.meta; + end + H5G.close(rootgid); + end + else + [varargout{1:nargout}] = load_one(loc, opt); + end + H5F.close(loc); +catch ME + H5F.close(loc); + rethrow(ME); +end + +if (jsonopt('Regroup', 0, opt)) + if (nargout >= 1) + varargout{1} = regrouph5(varargout{1}); + elseif (nargout >= 2) + varargout{2} = regrouph5(varargout{2}); + end +end + +if (isfield(opt, 'jdata') && opt.jdata && nargout >= 1) + varargout{1} = jdatadecode(varargout{1}, 'Base64', 0, opt); +end + +% -------------------------------------------------------------------------- +function [data, meta] = load_one(loc, opt) + +data = struct(); +meta = struct(); +inputdata = struct('data', data, 'meta', meta, 'opt', opt); + +% Load groups and datasets +try + [status, count, inputdata] = H5L.iterate(loc, opt.order, 'H5_ITER_INC', 0, @group_iterate, inputdata); +catch ME + if (strcmp(opt.order, 'H5_INDEX_CRT_ORDER')) + [status, count, inputdata] = H5L.iterate(loc, 'H5_INDEX_NAME', 'H5_ITER_INC', 0, @group_iterate, inputdata); + else + rethrow(ME); + end +end + +data = inputdata.data; +meta = inputdata.meta; + +% -------------------------------------------------------------------------- +function [status, res] = group_iterate(group_id, objname, inputdata) +status = 0; +attr = struct(); + +encodename = jsonopt('PackHex', 1, inputdata.opt); + +try + data = inputdata.data; + meta = inputdata.meta; + + % objtype index + info = H5G.get_objinfo(group_id, objname, 0); + objtype = info.type; + objtype = objtype + 1; + + if objtype == 1 + % Group + name = regexprep(objname, '.*/', ''); + + group_loc = H5G.open(group_id, name); + try + [sub_data, sub_meta] = load_one(group_loc, inputdata.opt); + H5G.close(group_loc); + catch ME + H5G.close(group_loc); + rethrow(ME); + end + if (encodename) + name = encodevarname(name); + else + name = genvarname(name); + end + data.(name) = sub_data; + meta.(name) = sub_meta; + + elseif objtype == 2 + % Dataset + name = regexprep(objname, '.*/', ''); + + dataset_loc = H5D.open(group_id, name); + try + sub_data = H5D.read(dataset_loc, ... + 'H5ML_DEFAULT', 'H5S_ALL', 'H5S_ALL', 'H5P_DEFAULT'); + try + [status, count, attr] = H5A.iterate(dataset_loc, 'H5_INDEX_NAME', 'H5_ITER_INC', 0, @getattribute, attr); + catch + attr = []; + end + H5D.close(dataset_loc); + catch exc + H5D.close(dataset_loc); + rethrow(exc); + end + if (ischar(sub_data) && numel(sub_data) > 1 && sub_data(end) == 0) + sub_data = sub_data(1:end - 1); + end + + if ((isnumeric(sub_data) && inputdata.opt.dotranspose) || (iscell(sub_data) && length(sub_data) > 1)) + sub_data = permute(sub_data, ndims(sub_data):-1:1); + end + sub_data = fix_data(sub_data, attr, inputdata.opt); + if (encodename) + name = encodevarname(name); + else + name = genvarname(name); + end + data.(name) = sub_data; + meta.(name) = attr; + end +catch ME + rethrow(ME); +end +res = struct('data', data, 'meta', meta, 'opt', inputdata.opt); + +% -------------------------------------------------------------------------- +function data = fix_data(data, attr, opt) +% Fix some common types of data to more friendly form. + +if isstruct(data) + fields = fieldnames(data); + + if (length(intersect(fields, {'SparseIndex', opt.complexformat{1}})) == 2) + if isnumeric(data.SparseIndex) && isnumeric(data.(opt.complexformat{1})) + if (nargin > 1 && isstruct(attr)) + if (isfield(attr, 'SparseArraySize')) + spd = sparse(1, prod(attr.SparseArraySize)); + if (isfield(data, opt.complexformat{2})) + spd(data.SparseIndex) = complex(data.(opt.complexformat{1}), data.(opt.complexformat{2})); + else + spd(data.SparseIndex) = data.(opt.complexformat{1}); + end + data = reshape(spd, attr.SparseArraySize(:)'); + return + end + end + end + else + if (numel(opt.complexformat) == 2 && length(intersect(fields, opt.complexformat)) == 2) + if isnumeric(data.(opt.complexformat{1})) && isnumeric(data.(opt.complexformat{2})) + data = data.(opt.complexformat{1}) + 1j * data.(opt.complexformat{2}); + end + else + % if complexformat is not specified or not found, try some common complex number storage formats + if (length(intersect(fields, {'Real', 'Imag'})) == 2) + if isnumeric(data.Real) && isnumeric(data.Imag) + data = data.Real + 1j * data.Imag; + end + elseif (length(intersect(fields, {'real', 'imag'})) == 2) + if isnumeric(data.real) && isnumeric(data.imag) + data = data.real + 1j * data.imag; + end + elseif (length(intersect(fields, {'Re', 'Im'})) == 2) + if isnumeric(data.Re) && isnumeric(data.Im) + data = data.Re + 1j * data.Im; + end + elseif (length(intersect(fields, {'re', 'im'})) == 2) + if isnumeric(data.re) && isnumeric(data.im) + data = data.re + 1j * data.im; + end + elseif (length(intersect(fields, {'r', 'i'})) == 2) + if isnumeric(data.r) && isnumeric(data.i) + data = data.r + 1j * data.i; + end + end + end + end +end + +if (isa(data, 'uint8') || isa(data, 'int8')) + if (nargin > 1 && isstruct(attr)) + if (isfield(attr, 'MATLABObjectClass')) + data = getArrayFromByteStream(data); % use undocumented function + end + end +end + +% handeling string arrays (or cell of char strings) +if (iscell(data) && length(data) > 1) + if (all(cellfun(@ischar, data)) && exist('string') && opt.stringarray) + data = string(data); + end +end + +% -------------------------------------------------------------------------- +function [status, dataout] = getattribute(loc_id, attr_name, info, datain) +status = 0; +attr_id = H5A.open(loc_id, attr_name, 'H5P_DEFAULT'); +datain.(attr_name) = H5A.read(attr_id, 'H5ML_DEFAULT'); +H5A.close(attr_id); +dataout = datain; diff --git a/loadjnifti.m b/loadjnifti.m new file mode 100644 index 0000000..84f5471 --- /dev/null +++ b/loadjnifti.m @@ -0,0 +1,57 @@ +function jnii = loadjnifti(filename, varargin) +% +% jnii=loadjnifti(inputfile) +% or +% jnii=loadjnifti(inputfile, 'Param1',value1, 'Param2',value2,...) +% +% Load a standard NIFTI-1/2 file or text or binary JNIfTI file with +% format defined in JNIfTI specification: https://github.com/NeuroJSON/jnifti +% +% author: Qianqian Fang (q.fang neu.edu) +% +% input: +% inputfile: the output file name to the JNIfTI or NIFTI-1/2 file +% *.bnii for binary JNIfTI file +% *.jnii for text JNIfTI file +% *.nii for NIFTI-1/2 files +% options: (optional) if loading from a .bnii file, please see the options for +% loadbj.m (part of JSONLab); if loading from a .jnii, please see the +% supported options for loadjson.m (part of JSONLab). +% +% output: +% jnii: a structure (array) or cell (array). The data structure can +% be completely generic or auxilary data without any JNIfTI +% constructs. However, if a JNIfTI object is included, it shall +% contain the below subfields (can appear within any depth of the +% structure) +% jnii.NIFTIHeader - a structure containing the 1-to-1 mapped NIFTI-1/2 header +% jnii.NIFTIData - the main image data array +% jnii.NIFTIExtension - a cell array contaiing the extension data buffers +% +% example: +% jnii=jnifticreate(uint8(magic(10)),'Name','10x10 magic matrix'); +% savejnifti(jnii, 'magic10.jnii') +% newjnii=loadjnifti('magic10.jnii'); +% +% this file is part of JNIfTI specification: https://github.com/NeuroJSON/jnifti +% +% License: Apache 2.0, see https://github.com/NeuroJSON/jnifti for details +% + +if (nargin < 1) + error('you must provide data and output file name'); +end + +if (~exist('savejson', 'file')) + error('you must first install JSONLab from http://github.com/fangq/jsonlab/'); +end + +if (regexp(filename, '\.nii$')) + jnii = nii2jnii(filename, 'jnii'); +elseif (regexp(filename, '\.jnii$')) + jnii = loadjson(filename, varargin{:}); +elseif (regexp(filename, '\.bnii$')) + jnii = loadbj(filename, varargin{:}); +else + error('file suffix must be .jnii for text JNIfTI, .bnii for binary JNIfTI or .nii for NIFTI-1/2 files'); +end diff --git a/loadnifti.m b/loadnifti.m new file mode 100644 index 0000000..1d1b533 --- /dev/null +++ b/loadnifti.m @@ -0,0 +1,21 @@ +function varargout = loadnifti (varargin) +% +% jnii=loadnifti(filename) +% or +% nii=loadnifti(filename,option) +% +% Read a NIfTI-1/2 (*.nii/.nii.gz) or Analyze 7.5 (*.hdr/*.img/.hdr.gz/.img.gz) +% image file. +% +% author: Qianqian Fang (q.fang neu.edu) +% +% Please run `help nii2jnii` to see the input output outputs. +% This function is an alias to nii2jnii +% +% +% this file is part of JNIfTI specification: https://github.com/NeuroJSON/jnifti +% +% License: Apache 2.0, see https://github.com/NeuroJSON/jnifti for details +% + +[varargout{1:nargout}] = nii2jnii(varargin{:}); diff --git a/memmapstream.m b/memmapstream.m new file mode 100644 index 0000000..b4497bb --- /dev/null +++ b/memmapstream.m @@ -0,0 +1,81 @@ +function outstruct = memmapstream(bytes, format, varargin) +% +% outstruct=memmapstream(bytes, format) +% +% Map a byte-array (in char array or uint8/int8 array) into a structure +% using a dictionary (format is compatible with memmapfile in MATLAB) +% +% This function is compatible with both MATLAB and GNU Octave. +% +% author: Qianqian Fang (q.fang neu.edu) +% +% input: +% bytes: a char, int8 or uint8 vector or array +% format: a 3-column cell array in the format compatible with the +% 'Format' parameter of memmapfile in MATLAB. It has the +% following structure +% +% column 1: data type string, it can be one of the following +% 'int8','int16','int32','int64', +% 'uint8','uint16','uint32','uint64', +% 'single','double','logical' +% column 2: an integer vector denoting the size of the data +% column 3: a string denoting the fieldname in the output struct +% +% For example format={'int8',[1,8],'key'; 'float',[1,1],'value'} +% reads the first 8 bytes from 'bytes' as the first subfield +% 'key' and the following 4 bytes as the floating point 'value' +% subfield. +% +% output: +% outstruct: a structure containing the required field +% +% example: +% bytestream=['Andy' 5 'JT']; +% format={'uint8', [1,4], 'name', +% 'uint8', [1,1], 'age', +% 'uint8', [1,2], 'school'}; +% data=memmapstream(bytestream,format); +% +% this file is part of JNIfTI specification: https://github.com/NeuroJSON/jnifti +% +% License: Apache 2.0, see https://github.com/NeuroJSON/jnifti for details +% + +if (nargin < 2) + error('must provide bytes and format as inputs'); +end + +if (~ischar(bytes) && ~isa(bytes, 'int8') && ~isa(bytes, 'uint8') || isempty(bytes)) + error('first input, bytes, must be a char-array or uint8/int8 vector'); +end + +if (~iscell(format) || size(format, 2) < 3 || size(format, 1) == 0 || ~ischar(format{1, 1})) + error('second input, format, must be a 3-column cell array, in a format described by the memmapfile Format field.'); +end + +bytes = bytes(:)'; + +datatype = struct('int8', 1, 'int16', 2, 'int32', 4, 'int64', 8, 'uint8', 1, 'uint16', 2, 'uint32', 4, 'uint64', 8, 'single', 4, 'double', 8); + +opt = varargin2struct(varargin{:}); +opt.usemap = jsonopt('usemap', 0, opt) && exist('containers.Map'); + +if (opt.usemap) + outstruct = containers.Map(); +else + outstruct = struct(); +end +len = 1; +for i = 1:size(format, 1) + bytelen = datatype.(format{i, 1}) * prod(format{i, 2}); + if (opt.usemap) + outstruct(format{i, 3}) = reshape(typecast(uint8(bytes(len:bytelen + len - 1)), format{i, 1}), format{i, 2}); + else + outstruct.(format{i, 3}) = reshape(typecast(uint8(bytes(len:bytelen + len - 1)), format{i, 1}), format{i, 2}); + end + len = len + bytelen; + if (len > length(bytes)) + break + end +end diff --git a/nii2jnii.m b/nii2jnii.m new file mode 100644 index 0000000..9a4580c --- /dev/null +++ b/nii2jnii.m @@ -0,0 +1,306 @@ +function nii = nii2jnii(filename, format, varargin) +% +% nii=nii2jnii(niifile, format, options) +% or +% nii2jnii(niifile, jniifile, options) +% nii=nii2jnii(niifile) +% +% A fast and portable NIFTI-1/2 and Analyze7.5 file parser and converter +% to the text and binary JNIfTI formats defined in JNIfTI specification: +% https://github.com/NeuroJSON/jnifti +% +% This function is compatible with both MATLAB and GNU Octave. +% It accepts .nii, .nii.gz, .hdr/.img and .hdr.gz/.img.gz input files +% +% author: Qianqian Fang (q.fang neu.edu) +% +% input: +% fname: the file name to the .nii, .nii.gz, .hdr/.img or .hdr.gz/.img.gz file +% format:'nii' for reading the NIfTI-1/2/Analyze files; +% 'jnii' to convert the nii data into an in-memory JNIfTI structure. +% 'niiheader' return only the nii header without the image data +% +% if format is not listed above and nii2jnii is called without +% an output, format must be a string specifying the output JNIfTI +% file name - *.jnii for text-based JNIfTI, or *.bnii for the +% binary version +% options: (optional) if saving to a .bnii file, please see the options for +% savebj.m (part of JSONLab); if saving to .jnii, please see the +% supported options for savejson.m (part of JSONLab). +% +% output: +% if the output is a JNIfTI data structure, it has the following subfield: +% nii.NIFTIHeader - a structure containing the 1-to-1 mapped NIFTI-1/2 header +% nii.NIFTIData - the main image data array +% nii.NIFTIExtension - a cell array contaiing the extension data buffers +% +% when calling as nii=nii2jnii(file,'nii'), the output is a NIFTI object containing +% 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 (for NIFTI-1) or 540 (for NIFTI-2) +% 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 +% +% dependency: +% No external dependency if reading .nii/.hdr/.img files; +% +% To load gzipped input files (.nii.gz/.hdr.gz/.img.gz), one must +% install the ZMat Toolbox (http://github.com/NeuroJSON/zmat) and +% JSONLab Toolbox (http://github.com/fangq/jsonlab); +% +% To save files into the text/binary JNIfTI formatted files, one +% need to install JSONLab (http://github.com/fangq/jsonlab). +% +% this file was initially developed for the MCX project: https://github.com/fangq/mcx/blob/master/utils/mcxloadnii.m +% +% this file is part of JNIfTI specification: https://github.com/NeuroJSON/jnifti +% +% License: Apache 2.0, see https://github.com/NeuroJSON/jnifti for details +% + +hdrfile = filename; +isnii = -1; +if (regexp(filename, '(\.[Hh][Dd][Rr](\.[Gg][Zz])*$|\.[Ii][Mm][Gg](\.[Gg][Zz])*$)')) + isnii = 0; +elseif (regexp(filename, '\.[Nn][Ii][Ii](\.[Gg][Zz])*$')) + isnii = 1; +end + +if (isnii < 0) + error('file must be a NIfTI (.nii/.nii.gz) or Analyze 7.5 (.hdr/.img,.hdr.gz/.img.gz) data file'); +end + +if (regexp(filename, '\.[Ii][Mm][Gg](\.[Gg][Zz])*$')) + hdrfile = regexprep(filename, '\.[Ii][Mm][Gg](\.[Gg][Zz])*$', '.hdr$1'); +end + +niftiheader = niiformat('nifti1'); + +if (~isempty(regexp(hdrfile, '\.[Gg][Zz]$', 'once')) || (exist('OCTAVE_VERSION', 'builtin') ~= 0)) + finput = fopen(hdrfile, 'rb'); + input = fread(finput, inf, 'uint8=>uint8'); + fclose(finput); + + if (regexp(hdrfile, '\.[Gg][Zz]$')) + if (~exist('gzipdecode', 'file')) + error('To process zipped files, you must install gzipdecode.m from the JSONLab toolbox: http://github.com/fangq/jsonlab'); + end + gzdata = gzipdecode(input); + else + gzdata = input; + end + clear input; + nii.hdr = memmapstream(gzdata, niftiheader); +else + fileinfo = dir(hdrfile); + if (isempty(fileinfo)) + error('specified file does not exist'); + end + header = memmapfile(hdrfile, ... + 'Offset', 0, ... + 'Writable', false, ... + 'Format', niftiheader(1:end - (fileinfo.bytes < 352), :)); + + nii.hdr = header.Data(1); +end + +[os, maxelem, dataendian] = computer; + +if (nii.hdr.sizeof_hdr ~= 348 && nii.hdr.sizeof_hdr ~= 540) + nii.hdr.sizeof_hdr = swapbytes(nii.hdr.sizeof_hdr); +end + +if (nii.hdr.sizeof_hdr == 540) % NIFTI-2 format + niftiheader = niiformat('nifti2'); + if (exist('gzdata', 'var')) + nii.hdr = memmapstream(gzdata, niftiheader); + else + header = memmapfile(hdrfile, ... + 'Offset', 0, ... + 'Writable', false, ... + 'Format', niftiheader(1:end - (fileinfo.bytes < 352), :)); + + nii.hdr = header.Data(1); + end +end + +if (nii.hdr.dim(1) > 7) + names = fieldnames(nii.hdr); + for i = 1:length(names) + nii.hdr.(names{i}) = swapbytes(nii.hdr.(names{i})); + end + if (nii.hdr.sizeof_hdr > 540) + nii.hdr.sizeof_hdr = swapbytes(nii.hdr.sizeof_hdr); + end + if (dataendian == 'B') + dataendian = 'little'; + else + dataendian = 'big'; + end +end + +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) % + 'single' 1 % float (32 bits/voxel) % + 'single' 2 % complex (64 bits/voxel) % + 'double' 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) % + 'int64' 1 % long long (64 bits) % + 'uint64' 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); +nii.endian = dataendian; + +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 + +if (nargin > 1 && strcmp(format, 'niiheader')) + return +end + +if (regexp(filename, '\.[Hh][Dd][Rr](\.[Gg][Zz])*$')) + filename = regexprep(filename, '\.[Hh][Dd][Rr](\.[Gg][Zz])*$', '.img$1'); +end + +imgbytenum = prod(double(nii.hdr.dim(2:nii.hdr.dim(1) + 1))) * nii.voxelbyte; + +if (isnii == 0 && ~isempty(regexp(filename, '\.[Gg][Zz]$', 'once'))) + finput = fopen(filename, 'rb'); + input = fread(finput, inf, 'uint8=>uint8'); + fclose(finput); + gzdata = gzipdecode(input); + nii.img = typecast(gzdata(1:imgbytenum), nii.datatype); +else + if (~exist('gzdata', 'var')) + fid = fopen(filename, 'rb'); + if (isnii) + fseek(fid, nii.hdr.vox_offset, 'bof'); + end + nii.img = fread(fid, imgbytenum, [nii.datatype '=>' nii.datatype]); + fclose(fid); + else + nii.img = typecast(gzdata(nii.hdr.vox_offset + 1:nii.hdr.vox_offset + imgbytenum), nii.datatype); + end +end + +nii.img = reshape(nii.img, nii.hdr.dim(2:nii.hdr.dim(1) + 1)); + +if (nargin > 1 && strcmp(format, 'nii')) + return +end + +nii0 = nii; + +nii = niiheader2jnii(nii0); + +nii.NIFTIData = nii0.img; + +if (isfield(nii0.hdr, 'extension') && nii0.hdr.extension(1) > 0) + if (exist('gzdata', 'var')) + nii.NIFTIExtension = cell(1); + count = 1; + bufpos = nii0.hdr.sizeof_hdr + 4; + while (bufpos < nii0.hdr.vox_offset) + nii.NIFTIExtension{count}.Size = typecast(gzdata(bufpos + 1:bufpos + 4), 'int32') - 8; + nii.NIFTIExtension{count}.Type = typecast(gzdata(bufpos + 5:bufpos + 8), 'int32'); + bufpos = bufpos + 8; + if (strcmp(dataendian, 'big')) + nii.NIFTIExtension{count}.Size = swapbytes(nii.NIFTIExtension{count}.Size); + nii.NIFTIExtension{count}.Type = swapbytes(nii.NIFTIExtension{count}.Type); + end + if (bufpos + nii.NIFTIExtension{count}.Size <= nii0.hdr.vox_offset) + nii.NIFTIExtension{count}.x0x5F_ByteStream_ = gzdata(bufpos + 1:bufpos + nii.NIFTIExtension{count}.Size); + end + bufpos = bufpos + bufpos + nii.NIFTIExtension{count}.Size; + count = count + 1; + end + else + fid = fopen(filename, 'rb'); + fseek(fid, nii0.hdr.sizeof_hdr + 4, 'bof'); + nii.NIFTIExtension = cell(1); + count = 1; + while (ftell(fid) < nii0.hdr.vox_offset) + nii.NIFTIExtension{count}.Size = fread(fid, 1, 'int32=>int32') - 8; + nii.NIFTIExtension{count}.Type = fread(fid, 1, 'int32=>int32'); + if (strcmp(dataendian, 'big')) + nii.NIFTIExtension{count}.Size = swapbytes(nii.NIFTIExtension{count}.Size); + nii.NIFTIExtension{count}.Type = swapbytes(nii.NIFTIExtension{count}.Type); + end + if (ftell(fid) + nii.NIFTIExtension{count}.Size < nii0.hdr.vox_offset) + nii.NIFTIExtension{count}.x0x5F_ByteStream_ = fread(fid, nii.NIFTIExtension{count}.Size, 'uint8=>uint8'); + end + count = count + 1; + end + fclose(fid); + end +end + +if (nargout == 0 && strcmp(format, 'nii') == 0 && strcmp(format, 'jnii') == 0) + if (~exist('savejson', 'file')) + error('you must first install JSONLab from http://github.com/fangq/jsonlab/'); + end + if (regexp(format, '\.jnii$')) + savejson('', nii, 'FileName', format, varargin{:}); + elseif (regexp(format, '\.bnii$')) + savebj('', nii, 'FileName', format, varargin{:}); + else + error('file suffix must be .jnii for text JNIfTI or .bnii for binary JNIfTI'); + end +end diff --git a/niicodemap.m b/niicodemap.m new file mode 100644 index 0000000..a21d96b --- /dev/null +++ b/niicodemap.m @@ -0,0 +1,93 @@ +function newval = niicodemap(name, value) +% +% newval=niicodemap(name, value) +% +% Bi-directional conversion from NIFTI codes to human-readable JNIfTI +% header string values +% +% author: Qianqian Fang (q.fang neu.edu) +% +% input: +% name: a header name as a string, currently support the below nii +% headers: 'intent_code', 'slice_code', 'datatype', 'qform', +% 'sform' and 'xyzt_units' and their corresponding JNIfTI +% headers: +% 'Intent','SliceType','DataType','QForm','SForm','Unit' +% value:the current header value, if it is a code, newval will +% output the string version; if it is a string, newval will +% return the code +% +% output: +% newval: the converted header value +% +% For the detailed nii header codes, please see +% https://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1.h +% +% example: +% newval=niicodemap('slice_code', '') +% newval=niicodemap('datatype', 'uint64') +% newval=niicodemap('datatype', 2) +% +% this file was initially developed for the MCX project: https://github.com/fangq/mcx/blob/master/utils/mcxloadnii.m +% +% this file is part of JNIfTI specification: https://github.com/NeuroJSON/jnifti +% +% License: Apache 2.0, see https://github.com/NeuroJSON/jnifti for details +% + +% code to name look-up-table + +if (~exist('containers.Map')) + newval = value; + return +end + +lut.intent_code = containers.Map([0, 2:24 1001:1011 2001:2005], ... + {'', 'corr', 'ttest', 'ftest', 'zscore', 'chi2', 'beta', ... + 'binomial', 'gamma', 'poisson', 'normal', 'ncftest', ... + 'ncchi2', 'logistic', 'laplace', 'uniform', 'ncttest', ... + 'weibull', 'chi', 'invgauss', 'extval', 'pvalue', ... + 'logpvalue', 'log10pvalue', 'estimate', 'label', 'neuronames', ... + 'matrix', 'symmatrix', 'dispvec', 'vector', 'point', 'triangle', ... + 'quaternion', 'unitless', 'tseries', 'elem', 'rgb', 'rgba', 'shape'}); + +lut.slice_code = containers.Map(0:6, {'', 'seq+', 'seq-', 'alt+', 'alt-', 'alt2+', 'alt-'}); + +lut.datatype = containers.Map([0, 2, 4, 8, 16, 32, 64, 128, 256, 512, 768, 1024, 1280, 1536, 1792, 2048, 2304], ... + {'', 'uint8', 'int16', 'int32', 'single', 'complex64', 'double', 'rgb24', 'int8', ... + 'uint16', 'uint32', 'int64', 'uint64', 'double128', 'complex128', ... + 'complex256', 'rgba32' }); + +lut.xyzt_units = containers.Map([0:3 8 16 24 32 40 48], ... + {'', 'm', 'mm', 'um', 's', 'ms', 'us', 'hz', 'ppm', 'rad'}); + +lut.qform = containers.Map(0:4, {'', 'scanner', 'aligned', 'talairach', 'mni'}); + +lut.unit = lut.xyzt_units; +lut.sform = lut.qform; +lut.slicetype = lut.slice_code; +lut.intent = lut.intent_code; + +% inverse look up table + +tul.intent_code = containers.Map(values(lut.intent_code), keys(lut.intent_code)); +tul.slice_code = containers.Map(values(lut.slice_code), keys(lut.slice_code)); +tul.datatype = containers.Map(values(lut.datatype), keys(lut.datatype)); +tul.xyzt_units = containers.Map(values(lut.xyzt_units), keys(lut.xyzt_units)); +tul.qform = containers.Map(values(lut.qform), keys(lut.qform)); + +tul.sform = tul.qform; +tul.slicetype = tul.slice_code; +tul.intent = tul.intent_code; +tul.unit = tul.xyzt_units; + +% map from code to name, or frmo name to code + +if (~isfield(lut, lower(name))) + error('property can not be found'); +end +if (~(ischar(value) || isa(value, 'string'))) + newval = lut.(lower(name))(value); +else + newval = tul.(lower(name))(value); +end diff --git a/niiformat.m b/niiformat.m new file mode 100644 index 0000000..20bd87c --- /dev/null +++ b/niiformat.m @@ -0,0 +1,128 @@ +function niiheader = niiformat(format) +% +% niiheader=niiformat(format) +% +% Return a NIfTI header format descriptor as an Nx3 cell array +% +% author: Qianqian Fang (q.fang neu.edu) +% +% input: +% format:'nifti1' - return the header descriptor for NIfTI-1 format +% 'nifti2' - return the header descriptor for NIfTI-2 format +% +% output: +% niiheader: an Nx3 cell array in the format similar to the 'Format' +% specifier of the memmapfile.m function in MATLAB +% It has the following structure: +% +% column 1: data type string, it can be one of the following +% 'int8','int16','int32','int64', +% 'uint8','uint16','uint32','uint64', +% 'single','double' +% column 2: an integer vector denoting the size of the data +% column 3: a string denoting the fieldname in the output struct +% +% this file is part of JNIfTI specification: https://github.com/NeuroJSON/jnifti +% +% License: Apache 2.0, see https://github.com/NeuroJSON/jnifti for details +% + +header.nifti1 = { ... + '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 % ... + }; + +header.nifti2 = { ... + 'int32' [1 1] 'sizeof_hdr' % !< MUST be 540 % % int sizeof_hdr; % ... + 'int8' [1 8] 'magic' % !< MUST be "ni2\0" or "n+2\0". % ... + 'int16' [1 1] 'datatype' % !< Defines data type! % % short datatype; % ... + 'int16' [1 1] 'bitpix' % !< Number bits/voxel. % % short bitpix; % ... + 'int64' [1 8] 'dim' % !< Data array dimensions.% % short dim[8]; % ... + 'double' [1 1] 'intent_p1' % !< 1st intent parameter. % % short unused8/9; % ... + 'double' [1 1] 'intent_p2' % !< 2nd intent parameter. % % short unused10/11; % ... + 'double' [1 1] 'intent_p3' % !< 3rd intent parameter. % % short unused12/13; % ... + 'double' [1 8] 'pixdim' % !< Grid spacings. % % float pixdim[8]; % ... + 'int64' [1 1] 'vox_offset' % !< Offset into .nii file % % float vox_offset; % ... + 'double' [1 1] 'scl_slope' % !< Data scaling: slope. % % float funused1; % ... + 'double' [1 1] 'scl_inter' % !< Data scaling: offset. % % float funused2; % ... + 'double' [1 1] 'cal_max' % !< Max display intensity % % float cal_max; % ... + 'double' [1 1] 'cal_min' % !< Min display intensity % % float cal_min; % ... + 'double' [1 1] 'slice_duration' % !< Time for 1 slice. % % float compressed; % ... + 'double' [1 1] 'toffset' % !< Time axis shift. % % float verified; % ... + 'int64' [1 1] 'slice_start' % !< First slice index. % % short dim_un0; % ... + 'int64' [1 1] 'slice_end' % !< Last slice index. % % float funused3; % ... + 'int8' [1 80] 'descrip' % !< any text you like. % % char descrip[80]; % ... + 'int8' [1 24] 'aux_file' % !< auxiliary filename. % % char aux_file[24]; % ... + 'int32' [1 1] 'qform_code' % !< NIFTI_XFORM_* code. % %-- all ANALYZE 7.5 --- % ... + 'int32' [1 1] 'sform_code' % !< NIFTI_XFORM_* code. % %below here are replaced% ... + 'double' [1 1] 'quatern_b' % !< Quaternion b param. % ... + 'double' [1 1] 'quatern_c' % !< Quaternion c param. % ... + 'double' [1 1] 'quatern_d' % !< Quaternion d param. % ... + 'double' [1 1] 'qoffset_x' % !< Quaternion x shift. % ... + 'double' [1 1] 'qoffset_y' % !< Quaternion y shift. % ... + 'double' [1 1] 'qoffset_z' % !< Quaternion z shift. % ... + 'double' [1 4] 'srow_x' % !< 1st row affine transform. % ... + 'double' [1 4] 'srow_y' % !< 2nd row affine transform. % ... + 'double' [1 4] 'srow_z' % !< 3rd row affine transform. % ... + 'int32' [1 1] 'slice_code' % !< Slice timing order. % ... + 'int32' [1 1] 'xyzt_units' % !< Units of pixdim[1..4] % ... + 'int32' [1 1] 'intent_code' % !< NIFTI_INTENT_* code. % % short unused14; % ... + 'int8' [1 16] 'intent_name' % !< 'name' or meaning of data. % ... + 'int8' [1 1] 'dim_info' % !< MRI slice ordering. % % char hkey_un0; % ... + 'int8' [1 15] 'reserved' % !< unused buffer % ... + 'int8' [1 4] 'extension' % !< header extension % ... + }; + +if (nargin < 1) + format = 'nifti1'; +end + +format = lower(format); + +if (isfield(header, format)) + niiheader = header.(format); +else + error('format must be either nifti1 or nifti2'); +end diff --git a/niiheader2jnii.m b/niiheader2jnii.m new file mode 100644 index 0000000..d5943c1 --- /dev/null +++ b/niiheader2jnii.m @@ -0,0 +1,65 @@ +function nii = niiheader2jnii(nii0) + +nii = struct(); +nii.NIFTIHeader.NIIHeaderSize = nii0.hdr.sizeof_hdr; +if (isfield(nii0.hdr, 'data_type')) + nii.NIFTIHeader.A75DataTypeName = deblank(char(nii0.hdr.data_type)); + nii.NIFTIHeader.A75DBName = deblank(char(nii0.hdr.db_name)); + nii.NIFTIHeader.A75Extends = nii0.hdr.extents; + nii.NIFTIHeader.A75SessionError = nii0.hdr.session_error; + nii.NIFTIHeader.A75Regular = nii0.hdr.regular; +end +nii.NIFTIHeader.DimInfo.Freq = bitand(uint8(nii0.hdr.dim_info), 7); +nii.NIFTIHeader.DimInfo.Phase = bitand(bitshift(uint8(nii0.hdr.dim_info), -3), 7); +nii.NIFTIHeader.DimInfo.Slice = bitand(bitshift(uint8(nii0.hdr.dim_info), -6), 7); +nii.NIFTIHeader.Dim = nii0.hdr.dim(2:2 + nii0.hdr.dim(1) - 1); +nii.NIFTIHeader.Param1 = nii0.hdr.intent_p1; +nii.NIFTIHeader.Param2 = nii0.hdr.intent_p2; +nii.NIFTIHeader.Param3 = nii0.hdr.intent_p3; +nii.NIFTIHeader.Intent = niicodemap('intent', nii0.hdr.intent_code); +nii.NIFTIHeader.DataType = niicodemap('datatype', nii0.hdr.datatype); +nii.NIFTIHeader.BitDepth = nii0.hdr.bitpix; +nii.NIFTIHeader.FirstSliceID = nii0.hdr.slice_start; +nii.NIFTIHeader.VoxelSize = nii0.hdr.pixdim(2:2 + nii0.hdr.dim(1) - 1); +nii.NIFTIHeader.Orientation = struct('x', 'r', 'y', 'a', 'z', 's'); +if (nii0.hdr.pixdim(1) < 0) + nii.NIFTIHeader.Orientation = struct('x', 'l', 'y', 'a', 'z', 's'); +end +nii.NIFTIHeader.NIIByteOffset = nii0.hdr.vox_offset; +nii.NIFTIHeader.ScaleSlope = nii0.hdr.scl_slope; +nii.NIFTIHeader.ScaleOffset = nii0.hdr.scl_inter; +nii.NIFTIHeader.LastSliceID = nii0.hdr.slice_end; +nii.NIFTIHeader.SliceType = niicodemap('slicetype', nii0.hdr.slice_code); +nii.NIFTIHeader.Unit.L = niicodemap('unit', bitand(uint8(nii0.hdr.xyzt_units), 7)); +nii.NIFTIHeader.Unit.T = niicodemap('unit', bitand(uint8(nii0.hdr.xyzt_units), 56)); +nii.NIFTIHeader.MaxIntensity = nii0.hdr.cal_max; +nii.NIFTIHeader.MinIntensity = nii0.hdr.cal_min; +nii.NIFTIHeader.SliceTime = nii0.hdr.slice_duration; +nii.NIFTIHeader.TimeOffset = nii0.hdr.toffset; +if (isfield(nii0.hdr, 'glmax')) + nii.NIFTIHeader.A75GlobalMax = nii0.hdr.glmax; + nii.NIFTIHeader.A75GlobalMin = nii0.hdr.glmin; +end +nii.NIFTIHeader.Description = deblank(char(nii0.hdr.descrip)); +nii.NIFTIHeader.AuxFile = deblank(char(nii0.hdr.aux_file)); +nii.NIFTIHeader.QForm = nii0.hdr.qform_code; +nii.NIFTIHeader.SForm = nii0.hdr.sform_code; +nii.NIFTIHeader.Quatern.b = nii0.hdr.quatern_b; +nii.NIFTIHeader.Quatern.c = nii0.hdr.quatern_c; +nii.NIFTIHeader.Quatern.d = nii0.hdr.quatern_d; +nii.NIFTIHeader.QuaternOffset.x = nii0.hdr.qoffset_x; +nii.NIFTIHeader.QuaternOffset.y = nii0.hdr.qoffset_y; +nii.NIFTIHeader.QuaternOffset.z = nii0.hdr.qoffset_z; +nii.NIFTIHeader.Affine(1, :) = nii0.hdr.srow_x; +nii.NIFTIHeader.Affine(2, :) = nii0.hdr.srow_y; +nii.NIFTIHeader.Affine(3, :) = nii0.hdr.srow_z; +nii.NIFTIHeader.Name = deblank(char(nii0.hdr.intent_name)); +nii.NIFTIHeader.NIIFormat = deblank(char(nii0.hdr.magic)); +if (isfield(nii0.hdr, 'extension')) + nii.NIFTIHeader.NIIExtender = nii0.hdr.extension; +end +nii.NIFTIHeader.NIIQfac_ = nii0.hdr.pixdim(1); +nii.NIFTIHeader.NIIEndian_ = nii0.endian; +if (isfield(nii0.hdr, 'reserved')) + nii.NIFTIHeader.NIIUnused_ = nii0.hdr.reserved; +end diff --git a/regrouph5.m b/regrouph5.m new file mode 100644 index 0000000..954a3f2 --- /dev/null +++ b/regrouph5.m @@ -0,0 +1,134 @@ +function data = regrouph5(root, varargin) +% +% data=regrouph5(root) +% or +% data=regrouph5(root,type) +% data=regrouph5(root,{'nameA','nameB',...}) +% +% Processing a loadh5 restored data and merge "indexed datasets", whose +% names start with an ASCII string followed by a contiguous integer +% sequence number starting from 1, into a cell array. For example, +% datasets {data.a1, data.a2, data.a3} will be merged into a cell/struct +% array data.a with 3 elements. +% +% A single subfield .name1 will be renamed as .name. Items with +% non-contigous numbering will not be grouped. If .name and +% .name1/.name2 co-exist in the input struct, no grouping will be done. +% +% The grouped subfield will appear at the position of the first +% pre-grouped item in the original input structure. +% +% author: Qianqian Fang (q.fang neu.edu) +% +% input: +% root: the raw input HDF5 data structure (loaded from loadh5.m) +% type: if type is set as a cell array of strings, it restrict the +% grouping only to the subset of field names in this list; +% if type is a string as 'snirf', it is the same as setting +% type as {'aux','data','nirs','stim','measurementList'}. +% +% output: +% data: a reorganized matlab structure. +% +% example: +% a=struct('a1',rand(5),'a2','string','a3',true,'d',2+3i,'e',{'test',[],1:5}); +% regrouph5(a) +% saveh5(a,'test.h5'); +% rawdata=loadh5('test.h5') +% data=regrouph5(rawdata) +% +% this file is part of EasyH5 Toolbox: https://github.com/NeuroJSON/easyh5 +% +% License: GPLv3 or 3-clause BSD license, see https://github.com/NeuroJSON/easyh5 for details +% + +if (nargin < 1) + help regrouph5; + return +end + +dict = {}; +if (~isempty(varargin)) + if (ischar(varargin{1}) && strcmpi(varargin{1}, 'snirf')) + dict = {'aux', 'data', 'nirs', 'stim', 'measurementList'}; + elseif (iscell(varargin{1})) + dict = varargin{1}; + end +end + +data = struct; +if (isstruct(root)) + data = repmat(struct, size(root)); + names = fieldnames(root); + newnames = struct(); + firstpos = struct(); + + for i = 1:length(names) + item = regexp(names{i}, '^(.*\D)(\d+)$', 'tokens'); + if (~isempty(item) && str2double(item{1}{2}) ~= 0 && ~isfield(root, item{1}{1})) + if (~isfield(newnames, item{1}{1})) + newnames.(item{1}{1}) = str2double(item{1}{2}); + else + newnames.(item{1}{1}) = [newnames.(item{1}{1}), str2double(item{1}{2})]; + end + if (~isfield(firstpos, item{1}{1})) + firstpos.(item{1}{1}) = length(fieldnames(data(1))); + end + else + for j = 1:length(root) + if (isstruct(root(j).(names{i}))) + data(j).(names{i}) = regrouph5(root(j).(names{i})); + else + data(j).(names{i}) = root(j).(names{i}); + end + end + end + end + + names = fieldnames(newnames); + if (~isempty(dict)) + names = intersect(names, dict); + end + + for i = length(names):-1:1 + len = length(newnames.(names{i})); + idx = newnames.(names{i}); + if ((min(idx) ~= 1 || max(idx) ~= len) && len ~= 1) + for j = 1:len + dataname = sprintf('%s%d', names{i}, idx(j)); + for k = 1:length(root) + if (isstruct(root(k).(dataname))) + data(k).(dataname) = regrouph5(root(k).(dataname)); + else + data(k).(dataname) = root(k).(dataname); + end + end + end + pos = firstpos.(names{i}); + len = length(fieldnames(data)); + data = orderfields(data, [1:pos, len, pos + 1:len - 1]); + continue + end + for j = 1:length(data) + data(j).(names{i}) = cell(1, len); + end + idx = sort(idx); + for j = 1:len + for k = 1:length(root) + obj = root(k).(sprintf('%s%d', names{i}, idx(j))); + if (isstruct(obj)) + data(k).(names{i}){j} = regrouph5(obj); + else + data(k).(names{i}){j} = obj; + end + end + end + pos = firstpos.(names{i}); + len = length(fieldnames(data)); + data = orderfields(data, [1:pos, len, pos + 1:len - 1]); + try + data.(names{i}) = cell2mat(data.(names{i})); + catch + end + end +end diff --git a/transposemat.m b/transposemat.m new file mode 100644 index 0000000..04593a0 --- /dev/null +++ b/transposemat.m @@ -0,0 +1,40 @@ +function data = transposemat(input) +% +% data=transposemat(input) +% +% Iterate over struct/cell and transpose 2D or higher-dimensional numerical +% array to match Octave loaded HDF5 array elements with loadh5 default setting +% +% author: Qianqian Fang (q.fang neu.edu) +% +% input: +% name: a matlab variable, can be a cell, struct, containers.Map, numeric array or strings +% +% output: +% newname: the restored original string +% +% example: +% a=struct('a', ones(2,3), 'b', 'a string', 'c', uint8(zeros(2,3,4))); +% b=transposemat(a) +% +% this file is part of EasyH5 Toolbox: https://github.com/NeuroJSON/easyh5 +% +% License: GPLv3 or 3-clause BSD license, see https://github.com/NeuroJSON/easyh5 for details +% + +if (isstruct(input)) + data = structfun(@transposemat, input, 'UniformOutput', false); +elseif (iscell(input)) + data = cellfun(@transposemat, input, 'UniformOutput', 'false'); +elseif (isa(input, 'containers.Map')) + allkeys = keys(input); + for i = 1:length(allkeys) + input(allkeys(i)) = transposemat(allkeys(i)); + end +elseif (isnumeric(input) && (ndims(input) > 2 || all(size(input) > 1))) + data = permute(input, ndims(input):-1:1); +elseif (ischar(input) && ndims(input) == 2 && size(input, 1) == 1 && size(input, 2) > 1 && input(end) == ' ') + data = input(1:end - 1); +else + data = input; +end