% Write NHDR/NRRD files
%
% headerInfo = nhdr_nrrd_write(nrrdFileName, headerInfo, bWriteData)
%
% Inputs:
% * nrrdFilename (char): path to NHDR/NRRD file(s), either a detached header,
% a detached header pointing to detached data files or a nrrd standalone
% depending on what is specified in headerInfo
% * headerInfo (struct): structure containing all the required and optional
% NRRD fields, such as produced by nhdr_nrrd_read.m, as well as the data.
% List of accepted structure field names with corresponding NRRD fields in
% parentheses:
% - byteskip (byte skip: <desc>): scalar integer. For now, can only
% accept -1 with raw encoding and 0 with all other encodings;
% - content (content: <desc>): string;
% - data (~): numerical array containing the raw data;
% - datafiles (data file: <desc>): cell array of strings containing the
% paths to detached data files relative to the header file. This can
% be a simple string if there is only one detached data file;
% - detached_header (~): Matlab structure containing all the fields
% required to write valid detached NRRD files containing the data.
% All detached NRRD data files are assumed to have the same header;
% - dimension (dimension: <desc>): scalar integer;
% - encoding (encoding: <desc>): string;
% - kinds (kinds: <desc>): cell array of strings;
% - lineskip (line skip: <desc>): scalar integer. For now, can only be
% set to zero unless raw encoding is used and byte skip is set to -1
% (which cancels the effect of line skip altogether);
% - measurementframe (measurement frame: <desc>): matrix of size
% [Ns, Ns], where Ns in the number of space dimensions;
% - sizes (sizes: <desc>): matrix of size [1, d] or [d, 1] where d is the
% dimension;
% - space (space: <desc>): string;
% - spacedimension (space dimension: <desc>): scalar integer;
% - spacedirections (space directions: <desc>): cell array of strings
% for code compatibility, e.g.
% {'(1,0,0)' '(0,2,0.1)' 'none' '(0.1,0.1,1.1)'} for 3
% world space dimensions;
% - spacedirections_matrix (~): matrix of size [Ns, Ns] where Ns is the
% number of space dimensions. Should match the string description
% contained in spacedirections, if specified;
% - spaceorigin (space origin: <desc>): matrix of size [1, Ns] or [Ns, 1]
% where Ns is the number os space dimensions;
% - spaceunits (space units: <desc>): cell array of strings;
% - type (type: <desc>): string.
%
% - bvalue (bvalue:=<value>): floating-point scalar (nominal b-value);
% - modality (modality:=<value>): string
% - gradients (DWMRI_gradient_%04d:=<value>): matrix of size [Ng, 3]
% where Ng is the number of gradients in a PGSE-type
% diffusion-weighted MRI protocol. Note that the b-value associated
% with each gradient [g_x, g_y, g_z] is computed as
% b = b_nom*(g_x^2 + g_y^2 + g_z^2), where b_nom is the
% nominal b-value, so the gradients must be scaled accordingly;
% * bWriteData (bool): set to true to write the data specified in
% headerInfo.data, set to false to just write the header file without the
% data
%
% Outputs:
% * headerInfo (struct): same structure as that passed as argument with
% additional (potentially mandatory) fields if inferrable from the data,
% such as 'dimension', 'sizes', etc.
%
% The function creates all the required NHDR/NRRD file(s).
%
% See nhdr_nrrd_read.m for more information.
%
%
% % Format definition:
% http://teem.sourceforge.net/nrrd/format.html
%
% Date: August 2018
% Author: Gaetan Rensonnet
function headerInfo = nhdr_nrrd_write(nrrdFileName, headerInfo, bWriteData)
[path_main_header, ~, fext] = fileparts(nrrdFileName);
if isfield(headerInfo, 'datafiles')
assert(strcmpi(fext, '.nhdr'), sprintf('Invalid filename extension ''%s''. A detached nrrd header (headerInfo structure contains ''datafiles'' field) should have extension ''.nhdr''.\n',fext));
else
assert(strcmpi(fext, '.nrrd'), sprintf('Invalid filename extension ''%s''. A standalone nrrd file (no ''datafiles'' field in headerInfo structure) should have extension ''.nrrd''.\n',fext));
end
header_fnames = fieldnames(headerInfo); % cell array of strings
fnames_parsed = false(1,length(header_fnames)); % flags field names from header successfully parsed
% --- Checking availability of required information before opening files ---
% Guess as much from the data as possible
% Check required 'sizes' information
if isfield(headerInfo, 'sizes')
if isfield(headerInfo, 'data')
assert(nrrd_size_check(size(headerInfo.data), headerInfo.sizes),...
sprintf('Sizes mismatch: [%s] in headerInfo structure is not compatible with the size of the data array [%s].\n',...
sprintf('%d ',headerInfo.sizes), sprintf('%d ',size(headerInfo.data)) ));
end
else
assert( isfield(headerInfo, 'data'), sprintf('Missing required field ''sizes'' in headerInfo structure, cannot be deduced from data because no data was provided.\n'));
headerInfo.sizes = size(headerInfo.data);
end
% Check required 'dimension' information
if isfield(headerInfo, 'dimension')
assert( isequal(headerInfo.dimension, length(headerInfo.sizes)), sprintf('Dimension mismatch: %d in headerInfo structure is not compatible with detected data size [%s].\n',...
headerInfo.dimension, sprintf('%d ', headerInfo.sizes)));
else
headerInfo.dimension = length(headerInfo.sizes);
end
% Check required 'type' information
if isfield(headerInfo, 'type')
matlabdatatype = nrrd_getMatlabDataType(headerInfo.type);
if isfield(headerInfo, 'data')
assert( isequal(class(headerInfo.data), matlabdatatype), ...
sprintf('Type mismatch: %s in headerInfo structure is not compatible with the Matlab type of the data array (%s).\nYou may want to look at nrrd_getMatlabDataType and at Matlab''s cast functions.',headerInfo.type, class(headerInfo.data)));
end
else
assert( isfield(headerInfo, 'data'), sprintf('Missing required field ''type'' in headerInfo structure, cannot be deduced from data because no data was provided.\n'));
headerInfo.type = get_nrrd_datatype(class(headerInfo.data));
end
% Check required 'encoding' information
if ~isfield(headerInfo, 'encoding')
headerInfo.encoding = 'raw';
end
% Set endianness to current platform's endianness
[~,~,endian] = computer();
if (isequal(endian, 'B'))
headerInfo.endian = 'big';
else
headerInfo.endian = 'little';
end
fidw = fopen( nrrdFileName, 'w');
if fidw == -1
error('nhdr_nrrd_writer: could not open file %s.\n',nrrdFileName);
end
% if an error is thrown while parsing, close open file before exiting
try
% --- Traditional header and required fields ---
fprintf( fidw, 'NRRD0005\n' );
fprintf( fidw, '# Complete NRRD file format specification at:\n' );
fprintf( fidw, '# http://teem.sourceforge.net/nrrd/format.html\n' );
if isfield(headerInfo,'content')
fprintf( fidw, 'content: %s\n', headerInfo.content); % always optional
fnames_parsed(strcmpi('content', header_fnames)) = true;
end
fprintf( fidw, 'type: %s\n', headerInfo.type );
fnames_parsed(strcmpi('type', header_fnames)) = true;
fprintf( fidw, 'encoding: %s\n', headerInfo.encoding);
fnames_parsed(strcmpi('encoding', header_fnames)) = true;
fprintf( fidw, 'dimension: %d\n', headerInfo.dimension );
fnames_parsed(strcmpi('dimension', header_fnames)) = true;
fprintf( fidw, 'sizes: ' );
for iI=1:length( headerInfo.sizes )
fprintf( fidw, '%d', headerInfo.sizes(iI) );
if ( iI~=length( headerInfo.sizes ) )
fprintf( fidw, ' ' );
end
end
fprintf( fidw, '\n' );
fnames_parsed(strcmpi('sizes', header_fnames)) = true;
fprintf( fidw, 'endian: %s\n', headerInfo.endian );
fnames_parsed(strcmpi('endian', header_fnames)) = true;
% ---- Defining orientation (always optional) -----
define_orientation = 0;
if isfield(headerInfo,'space')
fprintf(fidw, 'space: %s\n', headerInfo.space );
fnames_parsed(strcmpi('space', header_fnames)) = true;
define_orientation = 1;
num_sp_dim = nrrd_getSpaceDimensions(headerInfo.space);
elseif isfield (headerInfo,'spacedimension')
fprintf(fidw, 'space dimension: %d\n', headerInfo.spacedimension);
fnames_parsed(strcmpi('spacedimension', header_fnames)) = true;
define_orientation = 1;
num_sp_dim = headerInfo.spacedimension;
end
% Space and space dimension fields cannot be specified simultaneously
if isfield(headerInfo,'space') && isfield(headerInfo, 'spacedimension')
error('The always optional fields ''space'' and ''spacedimension'' defining orientation cannot be specified simultaneously in headerInfo structure');
end
if define_orientation
% So far the number of spatial dimensions should be 3 or 4
% (this is important for writing space directions, space origin,
% measurement frame)
assert(num_sp_dim>=3, sprintf('Number of space dimensions should be at least 3 in current nrrd format, detected %d instead in headerInfo structure.\n',num_sp_dim));
% check for necessary field 'space direction' when orientation defined
if ~isfield(headerInfo, 'spacedirections')
error('Missing field in headerInfo structure: spacedirections should always be specified if orientation is defined via ''space'' or ''space directions'' ');
end
assert(iscell(headerInfo.spacedirections), ...
['spacedirections field in headerInfo should be a cell array of ' ...
'strings for code compatibility, e.g. {''(1,0,0)'' ''(0,2,0.1)'' ''none'' ''(0.1,0.1,1.1)''}' ...
' for 3 world space dimensions.'])
% Check format and clean up when possible
for i = 1:length(headerInfo.spacedirections)
none_match = regexp(headerInfo.spacedirections{i}, 'none', 'match');
if ~isempty(none_match)
% none for non space axis
if ~strcmpi(headerInfo.spacedirections{i}, 'none')
fprintf(['nhdr_nrrd_write WARNING: detected %s instead of ' ...
'expected none in header.spacedirections{%d}.\n Wrote clean version to output file(s).\n'],...
headerInfo.spacedirections{i}, i)
end
headerInfo.spacedirections{i} = 'none';
else
% numerical vector for space (world) axis
btw_parentheses = headerInfo.spacedirections{i}(1) == '(' ...
&& headerInfo.spacedirections{i}(end) == ')';
if ~btw_parentheses
fprintf(['nhdr_nrrd_write WARNING: space vector in header.spacedirections{%d} ' ...
'should contain comma-separated values enclosed in single parentheses.\n ' ...
'At least one parenthesis missing here. Wrote clean version to output file(s).\n'],...
i)
end
headerInfo.spacedirections{i} = ['(' regexprep(headerInfo.spacedirections{i}, '[()]', ''), ')'];
end
end
spacedir_str = strjoin(headerInfo.spacedirections, ' ');
% Write to file
fprintf(fidw, 'space directions: %s\n', spacedir_str);
% Flag as parsed
fnames_parsed(strcmpi('spacedirections', header_fnames)) = true;
% Data defining space direction matrix
SD_data = strrep(spacedir_str, 'none', '');
SD_data = extractNumbersWithout(SD_data, {'(',')',',', '"', ''''}); % this fails if non numerical entries were not previously removed
assert(length(SD_data) == num_sp_dim^2, sprintf('expected spacedirections field to contain %d numbers (square of %d, the world space dimension). Found %d instead.\nspacedirections should be a cell array of strings containing vectors of the form (dx,dy,dz) or none entries.\n',...
num_sp_dim^2, num_sp_dim, length(SD_data)))
% Check spacedirections_matrix field (internal, not a NRRD field)
if isfield(headerInfo, 'spacedirections_matrix')
assert(isequal(headerInfo.spacedirections_matrix(:), SD_data(:)), ...
sprintf('Numeric data in spacedirections and spacedirections_matrix fields of headerInfo do not match. They should contain %d identical numbers in the same order.\n', num_sp_dim^2))
else
% add it to output headerInfo structure
headerInfo.spacedirections_matrix = reshape(SD_data(:), [num_sp_dim, num_sp_dim]);
end
fnames_parsed(strcmpi('spacedirections_matrix', header_fnames)) = true;
end
% --- optional options for the optional definition of orientation ---
% space origin
if isfield(headerInfo,'spaceorigin')
if define_orientation
so = headerInfo.spaceorigin;
assert(length(so)==num_sp_dim, ...
sprintf('Field ''spaceorigin'' in headerInfo structure expected vector with %d entries to match the defined orientation, detected %d entries instead.\n',num_sp_dim, length(so)));
fprintf( fidw, 'space origin: (%f%s)\n', so(1), sprintf(',%f',so(2:end)) );
fnames_parsed(strcmpi('spaceorigin', header_fnames)) = true;
else
error('Field ''spaceorigin'' in headerInfo structure cannot be specified if neither ''space'' nor ''spacedimension'' was specified.');
end
end
% measurement frame
if isfield(headerInfo,'measurementframe')
if define_orientation
mf = headerInfo.measurementframe;
assert(size(mf,1)==num_sp_dim && size(mf,2)==num_sp_dim,...
sprintf('Field ''measurementframe'' in headerInfo structure expected a %d-by-%d matrix to match the defined orientation, detected size [%s] instead.\n',...
num_sp_dim, num_sp_dim, sprintf('%d ',size(mf))));
fprintf( fidw, 'measurement frame:');
for imf = 1:size(mf,2)
fprintf(fidw,' (%f%s)', mf(1,imf), sprintf(',%f',mf(2:end,imf)));
end
fprintf(fidw, '\n');
fnames_parsed(strcmpi('measurementframe', header_fnames)) = true;
else
error('Field ''measurementframe'' in headerInfo structure cannot be specified if neither ''space'' nor ''spacedimension'' were specified.');
end
end
% spaceunits
if isfield(headerInfo,'spaceunits')
if define_orientation
fprintf( fidw, 'space units: ' );
for iI=1:length( headerInfo.spaceunits )
fprintf( fidw, '\"%s\"', char(headerInfo.spaceunits(iI)) );
if ( iI~=length( headerInfo.spaceunits ) )
fprintf( fidw, ' ' );
end
end
fprintf( fidw, '\n' );
fnames_parsed(strcmpi('spaceunits', header_fnames)) = true;
else
error('Field ''spaceunits'' in headerInfo structure cannot be specified if neither ''space'' nor ''spacedimension'' were specified.');
end
end
% --- end of orientation definition -----
% --- Optional fields "<field>: <desc>" ----
% kinds
if isfield(headerInfo, 'kinds')
assert(length(headerInfo.kinds)==headerInfo.dimension, sprintf('Detected %d kinds instead of expected nrrd dimension %d.\n',length(headerInfo.kinds), headerInfo.dimension));
fprintf( fidw, 'kinds: ' );
for iI=1:length( headerInfo.kinds )
fprintf( fidw, '%s', headerInfo.kinds{iI} );
if ( iI~=length( headerInfo.kinds ) )
fprintf( fidw, ' ' );
end
end
fprintf( fidw, '\n' );
fnames_parsed(strcmpi('kinds', header_fnames)) = true;
end
% line skip
if isfield(headerInfo, 'lineskip')
assert(headerInfo.lineskip >= 0, sprintf('Field ''lineskip'' in headerInfo structure should be non-negative, detected %d.\n', headerInfo.lineskip));
fprintf(fidw, 'lineskip: %d\n', headerInfo.lineskip);
fnames_parsed(strcmpi('lineskip', header_fnames)) = true;
end
% byte skip
if isfield(headerInfo, 'byteskip')
assert(headerInfo.byteskip >= -1, sprintf('Field ''byteskip'' should be -1 or a non-negative integer, detected %d.\n', headerInfo.byteskip));
if headerInfo.byteskip == -1
assert(strcmpi(headerInfo.encoding, 'raw'), sprintf('byte skip value of -1 is only valid with raw encoding.\n'));
end
fprintf(fidw, 'byteskip: %d\n', headerInfo.byteskip);
fnames_parsed(strcmpi('byteskip', header_fnames)) = true;
end
% --- "<key>:=<value>" pairs (always optional)
% Modality
if isfield(headerInfo, 'modality')
fprintf( fidw, 'modality:=%s\n',headerInfo.modality);
fnames_parsed( strcmpi('modality', header_fnames)) = true;
end
% b-value
if isfield(headerInfo, 'bvalue')
fprintf( fidw, 'DWMRI_b-value:=%9.8f\n', headerInfo.bvalue);
fnames_parsed( strcmpi('bvalue', header_fnames)) = true;
end
% DW-MRI gradient
if isfield(headerInfo, 'gradients')
Ngrads = size(headerInfo.gradients, 1);
for igrad = 1:Ngrads
strprint = strcat('DWMRI_gradient_',sprintf('%04d',igrad-1)) ; % The 0 flag in the %04d format specifier requests leading zeros in the output and sets minimum width of the printed value to 4
strprint = [strprint,':='] ;
strprint = [strprint, sprintf('%9.8f ', headerInfo.gradients(igrad,:))] ;
fprintf(fidw,[strprint '\n']) ;
end
fnames_parsed(strcmpi('gradients', header_fnames)) = true;
% FIXME: make it more general than numbering limited to 4 digits?
end
% --- External datafiles (should be performed in LIST mode) ---
if isfield(headerInfo, 'datafiles')
% This is a detached header, not a standalone nrrd file
% 'Convert' to cell array for code compatibility. This only works
% with one filename and fails with a comma-separated list of
% filenames for instance.
if ischar(headerInfo.datafiles)
headerInfo.datafiles = { headerInfo.datafiles};
end
if length(headerInfo.datafiles)==1
% data file: <filename>
fprintf(fidw, 'data file: %s\n',headerInfo.datafiles{1}); % path relative to detached header
else
% data file: LIST [<subdim>]
% FIXME: add support for subdim argument instead of ignoring it
fprintf(fidw, 'data file: LIST\n');
for i = 1:length(headerInfo.datafiles)
fprintf(fidw,'%s\n', headerInfo.datafiles{i});
end
end
fnames_parsed(strcmpi('datafiles', header_fnames)) = true;
end
catch me
fclose(fidw);
rethrow(me);
end
%% Write data
if bWriteData
if isfield(headerInfo, 'datafiles')
% This is a detached header, detached data files need to be
% written to
fclose(fidw); % close main file
N_data_tot = prod(headerInfo.sizes);
% 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;
for i = 1:N_data_files
% Check type of detached data file
[~,fname_data,ext_data] = fileparts(headerInfo.datafiles{i});
data_ind = (i-1)*N_data_per_file+1:i*N_data_per_file; % indices of data to be written
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
assert(isfield(headerInfo, 'detached_header'), sprintf('Missing field ''detached_header'' in headerInfo structure.\nThis is only required for detached headers with data stored in detached data files of type nrrd.\nShould contain a Matlab structure describing the chunk of data stored in the detached nrrd file (the data field, if present, will be overwritten).\nIt is assumed that the header is identical for all detached nrrd data files.\n'));
% we have to check for all i's because detached nrrd files may be
% interspersed with detached .raw or .hex files
info_detached_data = headerInfo.detached_header; % change headerInfo! put only part of the data, change dimensions, remove datalists
fnames_parsed(strcmpi('detached_header', header_fnames)) = true;
% Add data chunk to be written to new file
info_detached_data.data = headerInfo.data(data_ind);
info_detached_data.data = reshape(info_detached_data.data, info_detached_data.sizes); % for code compatibility
% Write detached nrrd file using a recursive call
bWrite_Detached_Data = true;
tmp_struct = nhdr_nrrd_write(fullfile(path_main_header, [fname_data, ext_data]), info_detached_data, bWrite_Detached_Data); % recursive call
else
% e.g., detached .raw or .hex file
fid_data = fopen( fullfile(path_main_header, [fname_data, ext_data]), 'w');
if( fid_data < 1 )
error('Detached data file number %d/%d (%s) could not be opened.\n', i, N_data_files, headerInfo.datafiles{i});
end
try
outtype = nrrd_getMatlabDataType(headerInfo.type);
writeData(fid_data, headerInfo.data(data_ind), outtype, headerInfo.encoding);
fclose(fid_data);
catch me_detached
fclose(fid_data);
rethrow(me_detached);
end
end
end
else
% Standalone nrrd file with data included after header
try
% After the header, there is a single blank line containing zero
% characters to separate it from data
fprintf(fidw, '\n');
outtype = nrrd_getMatlabDataType(headerInfo.type);
writeData(fidw, headerInfo.data, outtype, headerInfo.encoding);
fclose(fidw);
catch me_data
fclose(fidw);
rethrow(me_data);
end
end
else
fclose(fidw);
end
% Issue warnings for unknown/unsupported fields (after reading data)
for i = 1:length(header_fnames)
if ~fnames_parsed(i) && ~strcmpi(header_fnames{i},'data')
fprintf('nhdr_nrrd_write WARNING: could not parse and write field %s from headerInfo structure.\n', header_fnames{i});
end
end
end
% --- Auxiliary functions ------------
% ========================================================================
% 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
% ========================================================================
% Determine the nrrd datatype : from matlab datatype to outtype (NRRD)
% ========================================================================
function nrrddatatype = get_nrrd_datatype(matlab_metaType)
% Determine the datatype
switch (matlab_metaType)
case {'int8', 'uint8', 'int16', 'uint16', 'int32', 'uint32', 'int64',...
'uint64', 'double'}
nrrddatatype = matlab_metaType;
case {'single'}
nrrddatatype = 'float';
case {'block'}
error('Data type ''block'' is currently not supported.\n');
otherwise
error('Unknown matlab datatype %s', matlab_metaType)
end
end
% ========================================================================
% writeData -->
% fidIn is the open file we're overwriting
% matrix - data that have to be written
% datatype - type of data: int8, string, double...
% encoding - raw, gzip, txt/ascii
% ========================================================================
function ok = writeData(fidIn, matrix, datatype, encoding)
switch (encoding)
case {'raw'}
ok = fwrite(fidIn, matrix(:), datatype);
case {'gzip', 'gz'}
% Store in a raw file before compressing
tmpBase = tempname(pwd);
fidTmpRaw = fopen(tmpBase, 'w');
assert(fidTmpRaw > 3, 'Could not open temporary file for GZIP compression');
try
fwrite(fidTmpRaw, matrix(:), datatype);
fclose(fidTmpRaw);
catch me
fclose(fidTmpRaw);
delete(tmpBase);
rethrow(me);
end
% Now we gzip our raw file
tmpFile = [tmpBase '.gz'];
try
gzip(tmpBase); % this actually creates tmpFile
catch me
delete(tmpBase);
rethrow(me);
end
delete(tmpBase);
% Finally, we put this info into our nrrd file (fidIn)
fidTmpRaw = fopen(tmpFile, 'r'); % should this be in a try catch as well?
assert(fidTmpRaw > 3, 'Could not open temporary file for writing to nrrd file during gzip compression.');
try
% tmp = fread(fidTmpRaw, Inf, [datatype '=>' datatype]); % precision argument : from datatype (source) to datatype, how about just byte by byte ?
tmp = fread(fidTmpRaw, Inf, 'uint8=>uint8'); % this seems to be more robust
fclose(fidTmpRaw) ;
catch me
fclose(fidTmpRaw);
delete(tmpFile);
rethrow(me);
end
% ok = fwrite (fidIn, tmp, datatype); % why not byte by byte here?
ok = fwrite(fidIn, tmp, 'uint8'); % this seems to be more robust
delete (tmpFile);
case {'text', 'txt', 'ascii'}
ok = fprintf(fidIn,'%u ', matrix(:)); % FIX: better with %g ?
%ok = fprintf(fidIn,matrix(:), class(matrix));
otherwise
error('Unsupported encoding %s', encoding)
end
end
% ========================================================================
% Compare a Matlab-type size vector (no trailing ones, minimum length 2) to
% the size vector specified in a nrrd file/header which may contain
% trailing ones and may have only one dimension.
% ========================================================================
function sizesok = nrrd_size_check(matlab_size, nrrd_size)
% Make row vectors
matlab_size = matlab_size(:)';
nrrd_size = nrrd_size(:)';
% nrrd sizes may contain only one dimension
if length(nrrd_size)==1
sizesok = prod(matlab_size)==nrrd_size;
elseif length(nrrd_size)>=2
ind_max = find(nrrd_size~=1,1,'last'); % find last non-1 entry
% ones are ok in the first two dimensions:
if isempty(ind_max)
ind_max = 2;
else
ind_max = max(2, ind_max); % does not work if ind_max is empty
end
% Compare matlab size to nrrd size without trailing ones after from
% third entry on.
sizesok = isequal(matlab_size, nrrd_size(1:ind_max));
else
error('Argument nrrd_size is an empty matrix');
end
end