Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Matlab code to read in a NanoSIMS image file
function [im,pp,p]=read_im_file(ff,ask_for_planes,load_accumulated)
% Read the binary im file produced by the NanoSIMS 50L machine, and output
% the raw images as well as the additional parameters characterizing the
% dataset.
% NOTE: I don't know why, but data in some binary files are stored as uint8
% (e.g. from MPI Bremen), and in some others they are stored as ushort
% (e.g., IoW Warnemuende). It seems that this difference can be "detected"
% based on the value of the pixel_size in the Header_image structure.
% Therefore, the below coding relies on this to automatically find out how
% the binary data should be read.
% updated: L.P. 22-08-2012'
if nargin>1
afp=ask_for_planes;
else
afp=0;
end;
% some definitions
b8=2^8;
ttype = 'uint8';
num_bytes = 4;
% read the binary *.im file (header and the data afterwards)
disp(['Reading file ',ff,' ... ']);
fid=fopen(ff,'r');
% default flag for byte ordering, will be verified later
reverse_bytes = 0;
% default offset, will be modified later
offset = 452+8;
offset_old = offset;
% read the header size first, on byte 8. if the value is too large, which
% is indicative of reverse ordering of the bytes (i.e. that the file was
% generated under Windows), set the flag for the reverse byte ordering to
% 1. Also, change the offset value, which is used to read the masses in the
% binary file. (I dont have a clue why this offset depends on this, but it
% does, and I had to do it with a trial-and-error approach to find out the
% correct offset values for both types of data).
hsize = read_bytes(fid, 8, num_bytes, ttype, reverse_bytes, b8);
if hsize>2e6
reverse_bytes = ~reverse_bytes;
offset = 412+1*192+56;
end;
% read all the parameters in the header, at the positions specified i n the
% documentation
release = read_bytes(fid, 0, num_bytes, ttype, reverse_bytes, b8);
analysis_type = read_bytes(fid, 4, num_bytes, ttype, reverse_bytes, b8);
hsize = read_bytes(fid, 8, num_bytes, ttype, reverse_bytes, b8);
sample_type = read_bytes(fid, 12, num_bytes, ttype, reverse_bytes, b8);
data_present = read_bytes(fid, 16, num_bytes, ttype, reverse_bytes, b8);
sple_pos_x = read_bytes(fid, 20, num_bytes, ttype, reverse_bytes, b8);
sple_pos_y = read_bytes(fid, 24, num_bytes, ttype, reverse_bytes, b8);
analysis_name = read_bytes(fid, 28, 32, 'char', 0, b8);
username = read_bytes(fid, 60, 16, 'char', 0, b8);
samplename = read_bytes(fid, 76, 16, 'char', 0, b8);
ddate = remove_blanks(read_bytes(fid, 92, 16, 'char', 0, b8));
ttime = remove_blanks(read_bytes(fid, 108, 16, 'char', 0, b8));
filename = read_bytes(fid, 124, 16, 'char', 0, b8);
anal_duration = read_bytes(fid, 140, num_bytes, ttype, reverse_bytes, b8);
cycle_number = read_bytes(fid, 144, num_bytes, ttype, reverse_bytes, b8);
scantype = read_bytes(fid, 148, num_bytes, ttype, reverse_bytes, b8);
magnification = read_bytes(fid, 152, 2, ttype, reverse_bytes, b8);
sizetype = read_bytes(fid, 154, 2, ttype, reverse_bytes, b8);
size_detector = read_bytes(fid, 156, 2, ttype, reverse_bytes, b8);
no_used = read_bytes(fid, 158, 2, ttype, reverse_bytes, b8);
beam_blanking = read_bytes(fid, 160, num_bytes, ttype, reverse_bytes, b8);
pulverisation = read_bytes(fid, 164, num_bytes, ttype, reverse_bytes, b8);
pulve_duration = read_bytes(fid, 168, num_bytes, ttype, reverse_bytes, b8);
auto_cal_in_anal = read_bytes(fid, 172, num_bytes, ttype, reverse_bytes, b8);
autocal = read_bytes(fid, 176, 72, ttype, reverse_bytes, b8);
sig_reference = read_bytes(fid, 248, num_bytes, ttype, reverse_bytes, b8);
sigref = read_bytes(fid, 252, 156, ttype, reverse_bytes, b8);
nb_mass = read_bytes(fid, 408, num_bytes, ttype, reverse_bytes, b8);
% set to if 1 during debuging
if 0
ftmp=fopen([pwd delimiter 'tmp.txt'],'w');
fprintf(ftmp,'%f\n%f\n%f\n%f\n%f\n%f\n%f\n',release, analysis_type, hsize, ...
sample_type, data_present, sple_pos_x, sple_pos_y);
fprintf(ftmp,'%s\n%s\n%s\n%s\n%s\n%s\n',analysis_name, username, samplename, ...
ddate, ttime, filename);
fprintf(ftmp,'%d\n%d\n%d\n%d\n%d\n%d\n%d\n%d\n%d\n%d\n%d\n%d\n%d\n%d\n%d\n',...
anal_duration, cycle_number, scantype, magnification, sizetype, size_detector,...
no_used, beam_blanking, pulverisation, pulve_duration, auto_cal_in_anal, ...
autocal, sig_reference, sigref, nb_mass);
fclose(ftmp);
end;
% find the names and masses of the detected ions (I had to figure this out
% by trial and error, as the documentation did not really correspond to the
% final result).
at_mass = [];
mass_name = [];
for ii=0:(nb_mass-1)
[m1,h] = read_bytes(fid, offset+ii*192, 8, 'uint8', ~reverse_bytes, b8);
at_mass(1,ii+1) = typecast(uint8(h),'double');
[m1,h] = read_bytes(fid, offset+52+ii*192, 64, 'char', 0, b8);
if at_mass(1,ii+1)==0
mass_name{ii+1}='Esi';
else
if isempty(remove_blanks(m1))
% sometimes the element name is missing (such as in a file
% provided by Anand), so guess it or return AU (for atomic
% units)
switch round(at_mass(1,ii+1))
case 14, mn='14C';
case 33, mn='33S';
otherwise, mn=[num2str(round(at_mass(1,ii+1))) 'AU'];
end;
mass_name{ii+1} = mn;
else
mass_name{ii+1} = remove_blanks(m1);
end;
end;
end;
% for some unknown reason, the offset for reading the masses is not
% consistent between files produced in different laboratories, leading to
% at_mass(ii) = 0 for all ii. If this happens, read thd atomic masses
% again, with the old offset (reported by Kate Lewis)
if sum(at_mass)==0
fprintf(1,'*** Warning; Unusual input formatting. The originally detected masses may not be correct.\n');
fprintf(1,' Trying again with a different offset.\n');
for ii=0:(nb_mass-1)
[m1,h] = read_bytes(fid, offset_old+ii*192, 8, 'uint8', ~reverse_bytes, b8);
at_mass(1,ii+1) = typecast(uint8(h),'double');
[m1,h] = read_bytes(fid, offset_old+52+ii*192, 64, 'char', 0, b8);
if at_mass(1,ii+1)==0
mass_name{ii+1}='Esi';
else
mass_name{ii+1} = remove_blanks(m1);
end;
end;
end;
%
% read other properties of the images
width = read_bytes(fid, hsize-78, 2, 'uint8', reverse_bytes, b8);
height = read_bytes(fid, hsize-76, 2, 'uint8', reverse_bytes, b8);
pixel_size = read_bytes(fid, hsize-74, 2, 'uint8', reverse_bytes, b8);
n_images = read_bytes(fid, hsize-72, 2, 'uint8', reverse_bytes, b8);
n_planes = read_bytes(fid, hsize-70, 2, 'uint8', reverse_bytes, b8);
raster = read_bytes(fid, hsize-68, 4, 'uint8', reverse_bytes, b8);
nickname = read_bytes(fid, hsize-64, 64, 'char', 0, b8);
if ~load_accumulated
% this is used when the full binary data is going to be read
% this is where it is determined how to read the binary data
if pixel_size == 2
bin_format = 'uint8';
else
bin_format = 'ushort';
end;
% for some unknown reason, if the measurement was interrupted, n_planes may
% be equal to 0. in this case calculate the number of truly stored planes
% in the file from the number of masses and the size of the imaged area
if n_planes == 0
fileInfo = dir(ff);
n_planes=round((fileInfo.bytes-hsize)/(width*height*pixel_size*n_images));
end;
% true number of detected planes
Nptrue = n_planes;
fprintf(1,'Total number of detected planes: %d\n',Nptrue);
fprintf(1,'Size of the images: %d x %d pix\n',width,height);
fprintf(1,'Total number of detected masses: %d\n',length(mass_name));
% position of the 1st byte of the 1st plane of the 1st mass
p1 = hsize;
% ask for specific planes and masses that should be loaded from the im file
if afp
fprintf(1,'Detected masses:\n');
for i=1:length(mass_name)
fprintf(1,'%d:\t%s\n',i,mass_name{i});
end;
pm=input('Enter ID''s of the masses you want to load (e.g. [2 3 4], or empty for all): ');
pm_selected_empty=0;
if isempty(pm)
pm=[1:length(mass_name)];
pm_selected_empty=1;
else
ind = find(pm<=length(mass_name) & pm>=1);
pm = sort(pm(ind));
end;
b=[];
for i=1:length(pm)
b{i}=mass_name{pm(i)};
end;
mass_name = b;
at_mass = at_mass(pm);
pp=input('Enter planes you want to load (min:max, or empty for all): ');
if isempty(pp)
Np=Nptrue;
pp=[1:Np];
else
ind=find(pp<=Nptrue);
pp=pp(ind);
Np=length(pp);
end;
else
Np=Nptrue;
pp=[1:Np];
pm=1:nb_mass;
end;
fprintf(1,'Number of planes to be loaded: %d\n',Np);
% now read the image data for all selected masses and planes from byte p1
fseek(fid,p1,'bof');
kk=0;
N = width*height*2;
for ii=1:Nptrue
tmpim=[];
for jj=1:nb_mass
% read the image data in the specified bin_format (see above)
[d12,count]=fread(fid,N,bin_format);
% if ii==8
% a=1;
% end;
% get the higher and lower bytes and calculate the final image
d1=d12(1:2:N);
d2=d12(2:2:N);
if reverse_bytes
% when data have been stored under Windows
d=d2*b8+d1;
else
% when data have been stored under Unix
d=d1*b8+d2;
end;
a=vec2mat(d,width);
% for debugging, set this to 1 if you want to display every read image
if 0
figure(19);
imagesc(a);
set(gca,'dataaspectratio',[1 1 1]);
title(sprintf('plane=%d, mass=%d',ii,jj));
colormap(gray);
pause(0.5);
input('Press enter to continue');
end;
% fill the matrices of the masses
if nb_mass==1
tmpim = uint16(a); %double(a);
else
tmpim{jj} = uint16(a); %double(a);
end;
end;
% if the currently read plane is in the selected list of planes, and
% the currently mass in the selected list of masses,
% add it to im{jj}(:,:,:) at index kk, for all masses
if sum(ismember(pp,ii))>0
kk=kk+1;
%for jj=1:nb_mass
for jj=1:length(pm)
if nb_mass==1
im{jj}(:,:,kk) = tmpim;
else
im{jj}(:,:,kk) = tmpim{pm(jj)};
end;
end;
fprintf(1,'%d masses from plane %d (%d of %d) loaded.\n',length(pm),ii,kk,length(pp));
else
if ii==(pp(end)+1)
fprintf(1,'Please be patient until the im file is fully processed.\n');
end;
if ii==Nptrue
fprintf(1,'Done.\n');
end;
end;
% if ii<=pp(end) & kk>0
% if mod(kk,10)==0
% fprintf(1,'Plane %d (%d of %d) loaded.\n',pp(kk),kk,length(pp));
% elseif kk==length(pp)
% fprintf(1,'Plane %d (%d of %d) loaded.\n',pp(end),length(pp),length(pp));
% end;
% end;
end;
fclose(fid);
[pathstr, name, ext] = fileparts(ff);
if length(pp)==Nptrue & length(pm)==nb_mass
p.filename = [pathstr delimiter name];
else
s=[];
if length(pp)<Nptrue
s=sprintf('_p%d-%d',pp(1),pp(end));
end;
%if length(pm)<nb_mass
% s = [s '_m', sprintf('%d',pm)];
%end;
p.filename = sprintf('%s%c%s',pathstr,delimiter,name,s);
end;
else
% this is when only the accumulated data is going to be read (usually,
% this is useful when the processed data file is too huge to read
% everything at once)
fclose(fid);
[pathstr, name]=fileparts(ff);
indir = [pathstr delimiter name delimiter 'mat' delimiter];
loading_successful=1;
for i=1:nb_mass
infile=[indir mass_name{i} '.mat'];
if exist(infile,'file')
a=load(infile);
im_tmp{i}(:,:,1) = a.IM;
fprintf(1,'File %s loaded successfully.\n',infile);
else
loading_successful=0;
fprintf(1,'ERROR: File %s could not be loaded.\n',infile);
fprintf(1,'Please reprocess the RAW file to generate it.\n');
end;
end;
if loading_successful
im=im_tmp;
pp=1;
else
fprintf('ERROR: Some of the accumulated masses were not detected. Returning empty output.\n');
im=[];
pp=0;
end;
p.filename = [pathstr delimiter name];
end;
%fprintf(1,'Done.\n%d planes (%d-%d) stored.\n',kk,pp(1),pp(end));
% put the additional useful parameters to a structure p
p.reverse_bytes = reverse_bytes;
p.pos = [sple_pos_x sple_pos_y 0];
p.mass = mass_name;
p.mass_precise = at_mass;
p.date = ddate;
p.time = ttime;
p.width = width;
p.height = height;
p.scale = raster/1000;
a=0;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.