%INRLOAD loads an INR file
% This function allows the user to partially read a INR file in a matrix
% which is returned by the function.
% INRLOAD filename load the whole INR filename file.
% If additional arguments are specified, a specific part of the file is
% loaded:
% INRLOAD filename x_ini x_end only reads a column from all source images
% (from column x_ini to column x_end). First column is zero.
% INRLOAD filename x_ini x_end y_ini y_end only reads a square from
% all source images (from column x_ini to column x_end and from row x_ini
% to row x_end).
% INRLOAD filename x_ini x_end y_ini y_end z_ini z_end only reads a
% square from a specific interval of source images (from image z_ini to
% image z_end).
% INRLOAD filename x_ini x_end y_ini y_end z_ini z_end v_ini v_end
% only reads a square from a specific interval of source images of a
% specific interval of pixel dimensions (from dim. v_ini to dim. v_end).
% (in case a pixel is defined by a vector).
% Inf can be specified as coordinate value to designate the last
% coordinate in that dimension.
% im=INRLOAD returns a matrix im whose dimenions are size(im)=[X Y Z V].
% Take into account that when plotting images Matlab increases the first
% matrix dimension to plot a column. That is, the first matrix coordinate is
% the image height and the second one is the image width. Moreover, color
% channel dimension (V) is expected to be the third dimension. So, in order
% to show a grayscaled frame you first have to traspose it. Alternatively you can
% permute the dimensions of the whole INRLOAD output: permute(im,[2 1 4 3])
% [im,dl]=INRLOAD also returns a matrix dl of dimensions 2 x 4 in which
% each column represents the coordinate of the first and the last element
% loaded for each dimension.
%
% Usage example:
% a=inrload('sequence.inr', 0,Inf, 0,Inf, 0,99);
% it loads the first 100 complete frames from file sequence.inr.
%
% See also INRWRITE.
% Copyright (C) 2016 by Richard R. Carrillo
% $Revision: 1.3 $ $Date: 1/11/2016 $
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
function [voxels, dim_limits]=inrload(varargin)
nargs=nargin;
voxels=[]; % Defualt return values
dim_limits=[];
if nargs ~= 1 && nargs ~= 3 && nargs ~= 5 && nargs ~= 7 && nargs ~= 9
disp('You must specify 1, 3, 5, 7 or 9 arguments');
disp('inrload filename [x_ini x_end [y_ini y_end [z_ini z_end [v_ini v_end]]]]');
disp('All arguments must numeric except the first one')
disp('type: help inrload for more help');
else
filename = varargin{1};
fid = fopen(filename,'rb');
if fid ~= -1 % Successfully opened
header=inr_read_header(fid);
if isempty(header.read_error) % Header successfully read
% Display file data dimensions
disp(['Dimensions of data in file (X,Y,Z,V): ' num2str(header.n_dims(1)) sprintf(' x %i', header.n_dims(2:end))])
Dim_names='XYZV';
N_dims=length(Dim_names);
% If fn argument dimension data types are strings, convert to doubles
% The type may change denpending on how the fn is executed
for arg_ind=2:nargs
if isa(varargin{arg_ind},'char')
vargs_num(arg_ind-1) = str2double(varargin{arg_ind});
else
vargs_num(arg_ind-1) = varargin{arg_ind};
end
end
user_dim_size=[zeros(N_dims,1) ones(N_dims,1)*Inf]; % Default dimension coordinates of data to load: load all data: 0 -> Inf
% Replace dim_size values with the dimension limits specified by user
for arg_ind=1:(nargs-1)/2 % Parse input dimension coordenate intervals
user_dim_size(arg_ind,:) = [vargs_num(arg_ind*2 - 1) vargs_num(arg_ind*2)];
end
% If only 1 argument has been specified or [0,Inf] interval has
% been specified for all dimension, load the file completelly
if nargs == 1 || (all(user_dim_size(:,1)==0) && all(isinf(user_dim_size(:,2))))
disp(['loading file: ' filename ' completely']);
% Read the whole file at once to speed up the loading process
% Use the pixel data type specified in file header to read
% data from file and create voxels matrix with the same
% format (*)
voxels=fread(fid, Inf, ['*' header.matlab_data_type]);
disp(' 100%');
if numel(voxels)==prod(header.n_dims)
inr_dim_encoding=[header.n_dims(4) header.n_dims(1:3)]; % color hannels (V dim.) are encoded together
voxels=reshape(voxels, inr_dim_encoding); % convert vector obtained from fread into a matrix
voxels=permute(voxels,[2 3 4 1]); % Place the color channel dimension in the end as in CImg and INR header
%voxels=permute(voxels,[2 1 3 4]); % Permute X and Y dimensions to show images correctly with matlab image(voxels(:,:,1,1))
dim_limits = [zeros(1,N_dims) ; header.n_dims-1];
else
disp('Error: Number of voxel values in file does not correspond with dimension sizes in header')
end
else % only a part of the file must be loaded
dim_max_auto=isinf(user_dim_size(:,2)); % Dimensions in which the user want to load up to the last element
user_dim_size(dim_max_auto,2)=header.n_dims(dim_max_auto)-1;
dim_min_auto=isinf(user_dim_size(:,1)); % Dimensions in which the user want to load from the last element (that is, only the last element)
user_dim_size(dim_min_auto,1)=header.n_dims(dim_min_auto)-1;
stat_msg=['Partially loading activity file ' filename ': '];
for dim_ind=1:N_dims
stat_msg = [stat_msg Dim_names(dim_ind) '=[' num2str(user_dim_size(dim_ind,1)) ', ' num2str(user_dim_size(dim_ind,2)) '] '];
end
disp(stat_msg);
disp(' 00%')
% Check that input argument values are compatible with
% dimensions of data in the specified file
if all(user_dim_size >= 0)
if all(user_dim_size(:,1) < header.n_dims') && all(user_dim_size(:,2) < header.n_dims')
% Load binary (pixel) data from file
n_z_dim_end=(1+user_dim_size(3,2)-user_dim_size(3,1));
% Allocate matrix space (just for speed efficiency when adding values).
% Output matrix will has dimensions X,Y,Z,V (accoring to INR header field order and Cimg coordinates),
% that is voxels(width, height, frames, color_channels), however this is not the standard Matlab
% represetation for images/videos. Matlab expexts voxels(height, width, color_channels, frames),
% So, a matrix dimension permutation will be required after using this fn in order to use
% Matlab functions to show, play, etc the loaded data.
% But firstly we load frames as they are stored in file that is: V,X,Y,Z, so we allocate a matrix
% with these dimensions. Before exiting this fn we will permute dimensions to obtain CImg
% coordinates ordering in voxels matrix.
voxels=zeros(user_dim_size(4,2)-user_dim_size(4,1)+1, user_dim_size(1,2)-user_dim_size(1,1)+1, 1+user_dim_size(2,2)-user_dim_size(2,1), n_z_dim_end, header.matlab_data_type);
progress_end = n_z_dim_end; % For percentage display
process_update_period = ceil(progress_end / 100);
% Move file pointer to the first required data
% in Z dim, that is, to the first frame to
% load, which size is X*Y*V: prod(header.n_dims([1:2 4]))
fseek(fid, user_dim_size(3,1) * (header.data_size/8) * prod(header.n_dims([1:2 4])), 'cof');
for n_z_dim=1:n_z_dim_end
% For the remaining dimensions (X, Y and V) we load all the elements
% since it is normally faster than moving the file pointer many times
im_voxels=fread(fid, prod(header.n_dims([1:2 4])), ['*' header.matlab_data_type]); % Read a whole frame each time
% In an INR file color channels of the first
% pixel are stored first, then pixel of the
% first row and then the first frame, so
% dimensions: 4, 1, 2 and 3
im_voxels=reshape(im_voxels, header.n_dims([4 1:2]));
% Store the required part of the read frame in output matrix
voxels(:,:,:,n_z_dim) = im_voxels(1+(user_dim_size(4,1):user_dim_size(4,2)), 1+(user_dim_size(1,1):user_dim_size(1,2)), 1+(user_dim_size(2,1):user_dim_size(2,2)));
% Display progress percentage
current_progress = n_z_dim;
% we do not want to print the percentage so many times that we slow down the loading process
if mod(current_progress,process_update_period)==0 % Always true for progress_end < 100 since progress_end would be 1
fprintf(1,'\b\b\b\b% 3.f%%',100*current_progress/progress_end);
end
end
voxels=permute(voxels,[2 3 4 1]); % Set color channel as the last matrix dimension as in CImg and INR header
%voxels=permute(voxels,[2 1 3 4]); % Permute X and Y dimensions to show images correctly with matlab image(voxels(:,:,1,1))
dim_limits = user_dim_size'; % Return first and last coord. of loaded data
fprintf(1,'\b\b\b\b100%%\n');
else
disp('Error: all specified dimension limits must be lower than the dimensions in the specified file')
end
else
disp('Error: all specified dimension limits must be positive or zero')
end
end
if header.change_endianness
voxels=swapbytes(voxels); % File endianness and computer endianness are different, so data bytes must be swapped
end
else
disp(['Error reading file header: ' header.read_error])
end
fclose(fid);
else
disp(['Cannot open input file: ' varargin{1}]);
end
end
% Load INR file header and returns a struct containing the header fields
function header=inr_read_header(fid)
INR_header_ID='#INRIMAGE-4#';
INR_header_end='##}';
header=struct('read_error',''); % field read_ok will specify the obtained error or '' if the file header was successfully read
if fid ~= -1 % Valid file identifier
header_id=fscanf(fid,'%[^{]', 1); % Read format ID string
if strncmp(INR_header_ID,header_id,length(header_id)) == 1
fscanf(fid,' { \n', 1); % Skip block start and new line
end_header=false;
while ~end_header
header_line = fgetl(fid);
if ischar(header_line) % If a valid file line has been read:
if strcmp(INR_header_end, header_line) == 0 % If it is not the header end:
[param_name,param_name_read]=sscanf(header_line,'%[^=]=',1);
if param_name_read==1
param_value=sscanf(header_line,'%*[^=]=%[^\n]',1); % Skip param name and '=' and get param value
param_name=strtrim(param_name); % Remove leading and trailing space and tab character
param_value=strtrim(param_value);
switch param_name
case 'XDIM' % Length in X dimension
header.n_dims(1)=str2num(param_value);
case 'YDIM'
header.n_dims(2)=str2num(param_value);
case 'ZDIM'
header.n_dims(3)=str2num(param_value);
case 'VDIM'
header.n_dims(4)=str2num(param_value);
case 'VX'
header.voxel_size(1)=str2num(param_value);
case 'VY'
header.voxel_size(2)=str2num(param_value);
case 'VZ'
header.voxel_size(3)=str2num(param_value);
case 'TYPE'
header.data_type=param_value;
case 'PIXSIZE'
header.data_size=sscanf(param_value,' %u',1);
case 'SCALE'
header.scale=param_value;
case 'CPU'
header.cpu=param_value;
otherwise
if ~isempty(param_name) || ~isempty(param_value) % The line has a parameter name not recognized
warning(['Unexpected header parameter: ' param_name '=' param_value])
end
end
end
else
end_header=true; % header end found: exit loop
end
else
end_header=true; % file end found: exit
end
end
else
header.read_error='The format of specified file is unknown';
end
if isempty(header.read_error) % Header successfully read
% Translate INR data type name into Matlab data type name
if isfield(header, 'data_type') && isfield(header, 'data_size')
data_type_inr=[header.data_type '_' num2str(header.data_size)]; % Compose a data type name just to be used in switch statement
switch lower(data_type_inr)
case 'unsigned fixed_8'
header.matlab_data_type='uint8';
case 'signed fixed_8'
header.matlab_data_type='int8';
case 'unsigned fixed_16'
header.matlab_data_type='uint16';
case 'signed fixed_16'
header.matlab_data_type='int16';
case 'unsigned fixed_32'
header.matlab_data_type='uint32';
case 'signed fixed_32'
header.matlab_data_type='int32';
case 'unsigned fixed_64'
header.matlab_data_type='uint64';
case 'signed fixed_64'
header.matlab_data_type='int64';
case {'float_32','double_32'}
header.matlab_data_type='single';
case {'float_64','double_64'}
header.matlab_data_type='double';
otherwise
header.read_error='Incompatible data type specified in file header';
end
if isempty(header.read_error) % Data type successfully obtained
if isfield(header, 'cpu') % Endianness specified in header
[~,~,local_endianness]=computer;
switch header.cpu
case {'decm','alpha','pc'}
header.endianness='L';
case {'sun','sgi'}
header.endianness='B';
otherwise
header.read_error='Unknown CPU type specified in header';
end
if isfield(header, 'endianness')
header.change_endianness=(header.endianness~=local_endianness);
end
end
end
end
end
else
header.read_error='Invalid file identifier';
end