Skip to content
Snippets Groups Projects
Commit 853082aa authored by Nicola Vigano's avatar Nicola Vigano
Browse files

EDF: small changes 


Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent e4b3f91a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
zUtil_EDF/edf_check.m 0 → 100644
+ 18
− 0
View file @ 853082aa
function edf_check(filename)
% Acquire file information
edf_file = dir(filename);
% File Error Checking
if (isempty(edf_file))
error('EDF:file_not_found', 'Could not find: "%s"', filename);
end
if (edf_file.isdir == true)
error('EDF:file_is_directory', ...
'"%s" is a directory instead of a file.', filename);
end
if (edf_file.bytes == 0)
error('EDF:file_is_empty', '"%s" is an empty file.', filename);
end
end
\ No newline at end of file
zUtil_EDF/edf_info.m
+ 134
− 153
View file @ 853082aa
function edfheader = edf_info(fname)
function edf_header = edf_info(fname)
% EDF_INFO.M
% provides a structure containing all elements found in the header of the
% specified .edf file.
......@@ -12,7 +12,70 @@ function edfheader = edf_info(fname)
%
% Modified by Nicola Vigano, 07/12/2011, nicola.vigano@esrf.fr
function field_str = sfCleanField(field_str)
%%% Fields pattern
pattern = '(?<parameter>.*?)=(?<value>.*);';
edf_check(fname)
%%% Read the header
fid = fopen(fname, 'rt');
if (fid == -1)
error('EDF:file_not_found', ['Could not open ' fname]);
end
edf_header = [];
end_of_header = [];
edf_header.('headerlength') = 0;
while (isempty(end_of_header))
% Split lines
htxt = fgetl(fid);
% Check for inconsistencies in lines
if (isnumeric(htxt))
% Possible empty file! or at least malformed ;-)
% Close Resource
close(fid)
error('EDF:empty_line_in_header', 'Found an empty line in header.');
end
tmp = regexp(htxt, pattern, 'names');
if (~isempty(tmp))
field = sfCleanField(lower(strtrim(tmp.parameter)));
try
edf_header.(field) = strtrim(tmp.value);
catch mexc
gtPrintException(mexc)
warning('EDF:malformed_field', 'Didn''t understand this header');
disp(field);
end
end
% Look for end of header reached
end_of_header = find( htxt == '}' );
% Plus 1 is for the newline
edf_header.headerlength = edf_header.headerlength + length(htxt) +1;
end
% Count the bytes of the header, and check if it is a standard header
if (mod(edf_header.headerlength, 512) ~= 0)
warning('EDF:wrong_header_size', ...
'Header is not multiple of 512 bytes in %s (%d instead)!\n', ...
fname, edf_header.headerlength);
% disp('Setting correctly!')
% % XXX - may be cause of errors in edf_read. Check it!
% filePosition = round(filePosition/512) * 512;
end
% Release Resource
fclose(fid);
edf_header = sfCleanupHeaderInfo(edf_header);
end
function field_str = sfCleanField(field_str)
% if header field has weird characters (spaces, parentheses) - this removes
% them
field_str( find(field_str == ' ') ) = '_';
......@@ -20,162 +83,80 @@ function edfheader = edf_info(fname)
field_str( find(field_str == ',') ) = '_';
field_str( find(field_str == '(') ) = [];
field_str( find(field_str == ')') ) = [];
end
%% Fields pattern
pattern = '(?<parameter>.*?)=(?<value>.*);';
%% Acquire file information
edfFile = dir(fname);
%% File Error Checking
if isempty(edfFile)
Mexc = MException('EDF:fileNotFound', ['Could not find ' fname]);
throw(Mexc);
end
if edfFile.isdir == true
Mexc = MException('EDF:fileIsDirectory', [ fname ' is a directory instead of a file.']);
throw(Mexc);
end
if edfFile.bytes == 0
Mexc = MException('EDF:fileIsEmpty', [ fname ' is an empty file.']);
throw(Mexc);
end
%% Read the header
fid = fopen(fname, 'rt');
if fid == -1
Mexc = MException('EDF:fileNotFound', ['Could not open ' fname]);
throw(Mexc);
end
edfheader = [];
endofheader = [];
edfheader.('headerlength') = 0;
while isempty(endofheader)
% Split lines
htxt = fgetl(fid);
% Check for inconsistencies in lines
if isnumeric(htxt)
% Possible empty file! or at least malformed ;-)
% Close Resource
close(fid)
Mexc = MException('EDF:emptyLineInHeader', 'Found an empty line in header.');
throw(Mexc)
end
end
tmp = regexp(htxt, pattern, 'names');
if ~isempty(tmp)
field = sfCleanField(lower(strtrim(tmp.parameter)));
try
edfheader.(field) = strtrim(tmp.value);
catch Mexc
gtPrintException(Mexc)
warning('EDF:malformed_field', 'Didn''t understand this header');
disp(field);
end
function edf_header = sfCleanupHeaderInfo(edf_header)
% convert to numbers any fields that I know should not be strings
numeric_fields = { ...
'dim_1', 'dim_2', 'dim_3', 'size', 'count_time', 'timestamp_ms', ...
'col_beg', 'col_end', 'row_beg', 'row_end', 'col_bin', 'row_bin', ...
'energy', 'optic_used', 'acq_frame_nb', 'point_no', 'scan_no', ...
'run', 'image' ...
};
for ii = 1:length(numeric_fields)
if isfield(edf_header, numeric_fields{ii})
value = edf_header.(numeric_fields{ii});
value = str2double(value);
edf_header.(numeric_fields{ii}) = value;
end
end
% Look for end of header reached
endofheader = find( htxt == '}' );
% Plus 1 is for the newline
edfheader.headerlength = edfheader.headerlength + length(htxt) +1;
end
% Count the bytes of the header, and check if it is a standard header
if (mod(edfheader.headerlength, 512) ~= 0)
warning('EDF:wrong_header_size', ...
'Header is not multiple of 512 bytes in %s (%d instead)!\n', ...
fname, edfheader.headerlength);
% disp('Setting correctly!')
% % XXX - may be cause of errors in edf_read. Check it!
% filePosition = round(filePosition/512) * 512;
end
% Release Resource
fclose(fid);
edfheader = sfCleanupHeaderInfo(edfheader);
end
% convert all the motors and counters to individual fields
coupled_fields = { 'motor', 'counter' };
for ii = 1:length(coupled_fields)
field_pos = [coupled_fields{ii} '_pos'];
field_mne = [coupled_fields{ii} '_mne'];
if isfield(edf_header, field_pos)
field_tmp = [];
token_pos = textscan(edf_header.(field_pos), '%f');
token_name = textscan(edf_header.(field_mne), '%s');
for jj = 1:length(token_pos{1})
field_tmp.(token_name{1}{jj}) = token_pos{1}(jj);
end
edf_header.(coupled_fields{ii}) = field_tmp;
edf_header = rmfield(edf_header, {field_pos, field_mne} );
end
end
function edfheader = sfCleanupHeaderInfo(edfheader)
% convert to numbers any fields that I know should not be strings
numericFields = { ...
'dim_1', 'dim_2', 'dim_3', 'size', 'count_time', 'timestamp_ms', ...
'col_beg', 'col_end', 'row_beg', 'row_end', 'col_bin', 'row_bin', ...
'energy', 'optic_used', 'acq_frame_nb', 'point_no', 'scan_no', ...
'run', 'image' ...
};
for ii = 1:length(numericFields)
if isfield(edfheader, numericFields{ii})
value = edfheader.(numericFields{ii});
value = str2double(value);
edfheader.(numericFields{ii}) = value;
% convert endianness to standard representation
switch (edf_header.byteorder)
case 'HighByteFirst'
edf_header.byteorder = 'b';
case 'LowByteFirst'
edf_header.byteorder = 'l';
otherwise
disp('No endian-ness specified');
end
end
% convert all the motors and counters to individual fields
coupledFields = { 'motor', 'counter' };
for ii = 1:length(coupledFields)
field_pos = [coupledFields{ii} '_pos'];
field_mne = [coupledFields{ii} '_mne'];
if isfield(edfheader, field_pos)
field_tmp = [];
token_pos = textscan(edfheader.(field_pos), '%f');
token_name = textscan(edfheader.(field_mne), '%s');
for jj = 1:length(token_pos{1})
field_tmp.(token_name{1}{jj}) = token_pos{1}(jj);
end
edfheader.(coupledFields{ii}) = field_tmp;
edfheader = rmfield(edfheader, {field_pos, field_mne} );
% convert datatype to matlab standard representation
switch lower(edf_header.datatype)
case {'doublevalue'}
edf_header.datatype = 'float64';
case {'float','floatvalue','real'}
edf_header.datatype = 'float32';
case {'unsigned64'}
edf_header.datatype = 'uint64';
case {'signed64'}
edf_header.datatype = 'int64';
case {'unsignedinteger','unsignedlong'}
edf_header.datatype = 'uint32';
case {'signedinteger','signedlong'}
edf_header.datatype = 'int32';
case {'unsignedshort'}
edf_header.datatype = 'uint16';
case {'signedshort'}
edf_header.datatype = 'int16';
case {'unsignedbyte'}
edf_header.datatype = 'uint8';
case {'signedbyte'}
edf_header.datatype = 'int8';
otherwise
edf_header.datatype = [];
warning('EDF:unknown_type', ...
'Type "%s" not known', edf_header.datatype);
end
end
% convert endianness to standard representation
if strcmpi(edfheader.byteorder,'HighByteFirst')
edfheader.byteorder='b';
elseif strcmpi(edfheader.byteorder,'LowByteFirst')
edfheader.byteorder='l';
else
disp('No endian-ness specified');
end
% convert datatype to matlab standard representation
switch lower(edfheader.datatype)
case {'doublevalue'}
edfheader.datatype='float64';
case {'float','floatvalue','real'}
edfheader.datatype='float32';
case {'unsigned64'}
edfheader.datatype='uint64';
case {'signed64'}
edfheader.datatype='int64';
case {'unsignedinteger','unsignedlong'}
edfheader.datatype='uint32';
case {'signedinteger','signedlong'}
edfheader.datatype='int32';
case {'unsignedshort'}
edfheader.datatype='uint16';
case {'signedshort'}
edfheader.datatype='int16';
case {'unsignedbyte'}
edfheader.datatype='uint8';
case {'signedbyte'}
edfheader.datatype='int8';
otherwise
edfheader.datatype=[];
warning('EDF:unknown_type', ['Type ' edfheader.datatype ' not known']);
end
end
zUtil_EDF/edf_read.m
+ 132
− 153
View file @ 853082aa
function [img, varargout] = edf_read(filename, GT_bb, dispmes, info)
function [img, varargout] = edf_read(filename, GT_bb, dispmes, info, varargin)
% EDF_READ Reads images and 3D volumes from .edf files.
% --------------------------------------------------------------------------
% [img[, info]] = edf_read(filename[, GT_bb, dispmes, info]);
% [img[, info]] = edf_read(filename[, GT_bb, dispmes, info, varargin]);
%
% INPUT:
% filename = <string> Input edf filename
......@@ -20,170 +20,149 @@ function [img, varargout] = edf_read(filename, GT_bb, dispmes, info)
%
% Version 002 XX-06-2006 by NVigano, nicola.vigano@esrf.fr
% Many changes, many versions
% Modified by Nicola Vigano, 2011-2012, nicola.vigano@esrf.fr
% Modified by Nicola Vigano, 2011-2015, nicola.vigano@esrf.fr
%
% Version 001 XX-06-2006 by GJohnson
if (~exist('dispmes', 'var'))
dispmes = true;
else
dispmes = ~strcmp(dispmes, 'nodisp');
end
%% Acquire file information
edfFile = dir(filename);
%% File Error Checking
if (isempty(edfFile))
mexc = MException('EDF:fileNotFound', ['Could not find ' filename]);
throw(mexc);
end
if (edfFile.isdir == true)
mexc = MException('EDF:fileIsDirectory', [ filename ' is a directory instead of a file.']);
throw(mexc);
end
if (edfFile.bytes == 0)
mexc = MException('EDF:fileIsEmpty', [ filename ' is an empty file.']);
throw(mexc);
end
% Let's read the header
if (~exist('info', 'var') || isempty(info))
info = edf_info(filename);
elseif (~isfield(info, 'headerlength') || info.headerlength < 0)
% Make it sure the header length is found/computed!
info.headerlength = findHeaderLength(filename, info.byteorder);
end
% Test and fix third dimension
if (~isfield(info,'dim_3'))
info.dim_3 = 1;
else
if info.dim_3 > 1 && (exist('GT_bb','var') && ~isempty(GT_bb))
mexc = MException('EDF:roiIn3d', 'Cannot use an ROI on 3D EDF file');
throw(mexc);
if (~exist('dispmes', 'var'))
dispmes = true;
else
dispmes = ~strcmp(dispmes, 'nodisp');
end
end
fid = fopen(filename,'r',info.byteorder);
if (fid == -1)
mexc = MException('EDF:fileNotFound', ['Could not open ' fname]);
throw(mexc);
end
fseek(fid, info.headerlength, 'bof');
%%%%%%%%%%%%
%ROI
if exist('GT_bb','var') && ~isempty(GT_bb)
% GT_bb = [originx originy width height]
kodak_midline=1024-GT_bb(2)+1;
switch info.datatype
case {'uint8', 'int8'}
nbytes = 1;
case {'uint16', 'int16'}
nbytes = 2;
case {'uint32', 'int32', 'float32'}
nbytes = 4;
case {'uint64', 'int64', 'float64'}
nbytes = 8;
conf = struct('verbose', dispmes, 'convert', true);
conf = parse_pv_pairs(conf, varargin);
% Check file
edf_check(filename);
% Let's read the header
if (~exist('info', 'var') || isempty(info))
info = edf_info(filename);
elseif (~isfield(info, 'headerlength') || info.headerlength < 0)
% Make it sure the header length is found/computed!
info.headerlength = findHeaderLength(filename, info.byteorder);
end
% ROI was specified
pr_str = sprintf('%d*%s=>%s', GT_bb(3), info.datatype, info.datatype);
skip = nbytes * (info.dim_1 - GT_bb(3));
fseek(fid, nbytes*(info.dim_1*(GT_bb(2)-1)+GT_bb(1)-1), 0);
[img, count] = fread(fid, [GT_bb(3),GT_bb(4)], pr_str, skip);
dim_1 = GT_bb(3);
dim_2 = GT_bb(4);
else
%%%%%%%%%%%%%
kodak_midline = 1024;
dim_1 = info.dim_1;
dim_2 = info.dim_2;
[img, count] = fread(fid, info.dim_1*info.dim_2*info.dim_3, [info.datatype '=>' info.datatype]);
end
fclose(fid);
%% Check if the file was corrupted
if count ~= (dim_1*dim_2*info.dim_3)
report_string = [ 'File ' filename ' contains less data:\n' ...
' Expected: ' num2str(dim_1*dim_2*info.dim_3) ...
', Got: ' num2str(count)];
mexc = MException('EDF:fileCorrupted', report_string );
throw(mexc);
end
%% Let's reshape to the real format of the matrix
img = reshape(img, dim_1, dim_2, info.dim_3);
if (info.dim_3 > 1) % 3D volume
img = permute(img, [2 1 3]); % reorganise dimensions
else % 2D image
img = transpose(img);
end
if (~isa(img, 'double'))
img = double(img);
end
% Kodak4MV1: do any special processing required
if isfield(info,'ccdmodel') && any(isfield(info,{'isccdcorrected', 'israw'}))
% long test to keep compatible with 'israw' and 'isccdcorrected'
% fields. Should be able to remove 'israw' in due course
if strcmpi(info.ccdmodel,'kodak4mv1') ...
&& ( (isfield(info, 'isccdcorrected') && strcmpi(info.isccdcorrected,'false')) ...
|| (isfield(info, 'israw') && strcmpi(info.israw,'true')))
if kodak_midline > 0 && info.dim_2 == 2048
if dispmes == true
disp('Adding missing row for Kodak4M sensor and applying 0.96 gamma')
% Test and fix third dimension
if (~isfield(info, 'dim_3'))
info.dim_3 = 1;
else
if (info.dim_3 > 1 && (exist('GT_bb', 'var') && ~isempty(GT_bb)))
error('EDF:roi_in_3d', 'Cannot use an ROI on 3D EDF file');
end
imc = [ img(1:kodak_midline,:); ...
mean([img(kodak_midline,:); ...
img(kodak_midline+1,:)]); ...
img(kodak_midline+1:end-1,:) ];
img = imc .^ 0.96;
else
if (dispmes == true)
disp('WARNING: Cropped image may be corrupted due to Kodak4M sensor.')
end
end
end
fid = fopen(filename, 'r', info.byteorder);
if (fid == -1)
error('EDF:file_not_found', ['Could not open ' fname]);
end
end
%% If requested, let's output the header
if (nargout == 2)
varargout{1} = info;
end
function headerLength = findHeaderLength(filename, byteorder)
fidFind = fopen(filename, 'r', byteorder);
if (fidFind == -1)
mexcf = MException('EDF:fileNotFound', ['Could not open ' fname]);
throw(mexcf);
fseek(fid, info.headerlength, 'bof');
%%%%%%%%%%%%
%ROI
if (exist('GT_bb', 'var') && ~isempty(GT_bb))
% GT_bb = [originx originy width height]
kodak_midline = 1024 - GT_bb(2) + 1;
switch info.datatype
case {'uint8', 'int8'}
nbytes = 1;
case {'uint16', 'int16'}
nbytes = 2;
case {'uint32', 'int32', 'float32'}
nbytes = 4;
case {'uint64', 'int64', 'float64'}
nbytes = 8;
end
for ii = 1:8
fseek(fidFind, 512*ii -2, 'bof');
if (fgets(fidFind, 1) == '}')
fclose(fidFind);
headerLength = 512*ii;
return
% ROI was specified
pr_str = sprintf('%d*%s=>%s', GT_bb(3), info.datatype, info.datatype);
skip = nbytes * (info.dim_1 - GT_bb(3));
fseek(fid, nbytes * (info.dim_1 * (GT_bb(2) - 1) + GT_bb(1) - 1), 0);
[img, count] = fread(fid, [GT_bb(3), GT_bb(4)], pr_str, skip);
dim_1 = GT_bb(3);
dim_2 = GT_bb(4);
else
%%%%%%%%%%%%%
kodak_midline = 1024;
dim_1 = info.dim_1;
dim_2 = info.dim_2;
[img, count] = fread(fid, info.dim_1*info.dim_2*info.dim_3, [info.datatype '=>' info.datatype]);
end
fclose(fid);
%%% Check if the file was corrupted
if (count ~= (dim_1*dim_2*info.dim_3))
error('EDF:corrupted_file', ...
'File "%s" contains less data than expected:\n Expected: %d bytes, Got: %d bytes', ...
dim_1 * dim_2 * info.dim_3, count);
end
%%% Let's reshape to the real format of the matrix
img = reshape(img, dim_1, dim_2, info.dim_3);
if (info.dim_3 > 1) % 3D volume
img = permute(img, [2 1 3]); % reorganise dimensions
else % 2D image
img = transpose(img);
end
if (~isa(img, 'double') && conf.convert)
img = double(img);
end
% Kodak4MV1: do any special processing required
if (isfield(info, 'ccdmodel') && any(isfield(info, {'isccdcorrected', 'israw'})))
% long test to keep compatible with 'israw' and 'isccdcorrected'
% fields. Should be able to remove 'israw' in due course
if (strcmpi(info.ccdmodel, 'kodak4mv1') ...
&& ( (isfield(info, 'isccdcorrected') && strcmpi(info.isccdcorrected, 'false')) ...
|| (isfield(info, 'israw') && strcmpi(info.israw, 'true'))) )
if (kodak_midline > 0 && info.dim_2 == 2048)
if (conf.verbose)
disp('Adding missing row for Kodak4M sensor and applying 0.96 gamma')
end
imc = [ img(1:kodak_midline, :); ...
mean([img(kodak_midline, :); ...
img(kodak_midline+1, :)]); ...
img(kodak_midline+1:end-1, :) ];
img = imc .^ 0.96;
elseif (conf.verbose)
disp('WARNING: Cropped image may be corrupted due to Kodak4M sensor.')
end
end
end
%%% If requested, let's output the header
if (nargout == 2)
varargout{1} = info;
end
end
function header_length = findHeaderLength(filename, byteorder)
fid_find = fopen(filename, 'r', byteorder);
if (fid_find == -1)
error('EDF:file_not_found', ['Could not open ' fname]);
end
fprintf('Malformed header, searching manually.\n');
fseek(fidFind, 1, 'bof');
endofheader = false;
while (~endofheader)
htxt = fgetl(fidFind);
endofheader = find( htxt == '}' );
for ii = 1:8
fseek(fid_find, 512 * ii - 2, 'bof');
if (fgets(fid_find, 1) == '}')
fclose(fid_find);
header_length = 512 * ii;
return
end
headerLength = ftell(fidFind);
fclose(fidFind);
end
end % end of function
fprintf('Malformed header, searching manually.\n');
fseek(fid_find, 1, 'bof');
end_of_header = false;
while (~end_of_header)
htxt = fgetl(fid_find);
end_of_header = find( htxt == '}' );
end
header_length = ftell(fid_find);
fclose(fid_find);
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment