forked from flatironinstitute/NoRMCorre
-
Notifications
You must be signed in to change notification settings - Fork 0
/
bigread2.m
executable file
·207 lines (188 loc) · 8.03 KB
/
bigread2.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
function imData=bigread2(path_to_file,sframe,num2read)
%reads tiff files in Matlab bigger than 4GB, allows reading from sframe to sframe+num2read-1 frames of the tiff - in other words, you can read page 200-300 without rading in from page 1.
%based on a partial solution posted on Matlab Central (http://www.mathworks.com/matlabcentral/answers/108021-matlab-only-opens-first-frame-of-multi-page-tiff-stack)
%Darcy Peterka 2014, v1.0
%Darcy Peterka 2014, v1.1
%Darcy Peterka 2016, v1.2(bugs to [email protected])
%Eftychios Pnevmatikakis 2016, v1.3 (added hdf5 support)
%Program checks for bit depth, whether int or float, and byte order. Assumes uncompressed, non-negative (i.e. unsigned) data.
%
% Usage: my_data=bigread('path_to_data_file, start frame, num to read);
% "my_data" will be your [M,N,frames] array.
%Will do mild error checking on the inputs - last two inputs are optional -
%if nargin == 2, then assumes second is number of frames to read, and
%starts at frame 1
[~,~,ext] = fileparts(path_to_file);
if strcmpi(ext,'.tiff') || strcmpi(ext,'.tif');
%get image info
info = imfinfo(path_to_file);
% if ~isfield(info,'ImageDescription')
% blah=size(info);
% numFrames= blah(1);
% else
% he=info.ImageDescription;
% numFramesStr = regexp(he, 'images=(\d*)', 'tokens');
% numFrames = str2double(numFramesStr{1}{1});
% end
blah=size(info);%Nadav: Why not modulate this function by creating a small tool to count frames?
numFrames= blah(1);
num_tot_frames=numFrames;% why use another variable?
%should add more error checking for args... very ugly code below. works
%for me after midnight though...
if nargin<2
sframe = 1;
end
if nargin<3
num2read=numFrames-sframe+1;
end
if sframe<=0
sframe=1;
end
if num2read<1
num2read=1;
end
if sframe>num_tot_frames
sframe=num_tot_frames;
num2read=1;
display('starting frame has to be less than number of total frames...');
end
if (num2read+sframe<= num_tot_frames+1)
lastframe=num2read+sframe-1;
else
num2read=numFrames-sframe+1;
lastframe=num_tot_frames;
display('Hmmm...just reading from starting frame until the end');
end
bd=info.BitDepth;
he=info.ByteOrder;
bo=strcmp(he,'big-endian');
if (bd==64)
form='double';
elseif (bd==32)
form='single';
elseif (bd==16)
form='uint16';
elseif (bd==8)
form='uint8';
end
% Use low-level File I/O to read the file
fp = fopen(path_to_file , 'rb');
% The StripOffsets field provides the offset to the first strip. Based on
% the INFO for this file, each image consists of 1 strip.
he=info.StripOffsets;%nadav:seems to return only one number instead of an 1Xn array.
%Epnev:finds the offset of each strip in the movie. Image does not have to have
%uniform strips, but needs uniform bytes per strip/row.
idss=max(size(info(1).StripOffsets));%% - nadav: why is he taking the offset from the first image only?
%nadav:seems this function is ment to check the number of strips per frame, assuming that all frames has the same number of strips.
ofds=zeros(1,numFrames);%nadav:Creates an unnecesary 2000X2000 matrix, while in practice only a 2000X1 is requiered.%Nadav 07/10/17: fixed
for i=1:numFrames
ofds(i)=info(i).StripOffsets(1);
%ofds(i)
%nadav: This loop fills ofds matrix with the offset of the first strip in the frame.
%%nadav: It seems that the original reasoning for thisfor loop was to allow
%%for images with different number of stripes in each frame.
%%However, the variable 'idss' is computed with the assumption that the number
%%of strips per frame is uniform.
%%In our implementation, should we choose to assume that the number of stripes
%%per frame is uniform, the shorthand version ofds = [info.StripOffsets] would be faster
%%and more elegant.
end
sframemsg = ['Reading from frame ',num2str(sframe),' to frame ',num2str(num2read+sframe-1),' of ',num2str(num_tot_frames), ' total frames'];
disp(sframemsg)
pause(.2)
%go to start of first strip
fseek(fp, ofds(1), 'bof');
%framenum=numFrames;
framenum=num2read;
imData=cell(1,framenum);
he_w=info.Width;
he_h=info.Height;
% mul is set to > 1 for debugging only
mul=1;
if strcmpi(form,'uint16') || strcmpi(form,'uint8')
if bo
for cnt = sframe:lastframe
%cnt;
fseek(fp,ofds(cnt),'bof');%Nadav: Not sure why this function is needed since the file pointer moves tot the position after the file read.
tmp1 = fread(fp, [he_w he_h*mul], form, 0, 'ieee-be')';%nadav:why not use the fread argument 'skip' instead of fseek.
imData{cnt-sframe+1}=cast(tmp1,form);
end
else
for cnt=sframe:lastframe
% cnt;
fseek(fp,ofds(cnt),'bof');
tmp1 = fread(fp, [he_w he_h*mul], form, 0, 'ieee-le')';
imData{cnt-sframe+1}=cast(tmp1,form);
end
end
%Nadav: This code is very redundant. Suggested implementation is:
%for cnt = 1:lastframe-sframe
% tmp1 = fread(fp, [he_w he_h*mul], form, ofds(cnt), info.ByteOrder)';
%fileID , %dimensions of frame, %BitDepth, %offset, %ByteOrder
% imData{cnt}=cast(tmp1,form);
%end
%better implementation at the end of the if tree.
elseif strcmpi(form,'single')
if(bo)
for cnt = sframe:lastframe
%cnt;
fseek(fp,ofds(cnt),'bof');
tmp1 = fread(fp, [he_w he_h*mul], form, 0, 'ieee-be')';
imData{cnt-sframe+1}=cast(tmp1,'single');
end
else
for cnt = sframe:lastframe
%cnt;
fseek(fp,ofds(cnt),'bof');
tmp1 = fread(fp, [he_w he_h*mul], form, 0, 'ieee-le')';
imData{cnt-sframe+1}=cast(tmp1,'single');
end
end
elseif strcmpi(form,'double')
if(bo)
for cnt = sframe:lastframe
%cnt;
fseek(fp,ofds(cnt),'bof');
tmp1 = fread(fp, [he_w he_h*mul], form, 0, 'ieee-be.l64')';
imData{cnt-sframe+1}=cast(tmp1,'single');
end
else
for cnt = sframe:lastframe
%cnt;
fseek(fp,ofds(cnt),'bof');
tmp1 = fread(fp, [he_w he_h*mul], form, 0, 'ieee-le.l64')';
imData{cnt-sframe+1}=cast(tmp1,'single');
end
end
end
%This whole if tree is very redundant because a simple if statment setting form to single in case
%it is a double, would have solved this
%Suggested implementation:
%if strcmpi(form,'double')
% form = 'single'
% ByteOrder = ieee-be.l64;
%end
%for cnt = 1:lastframe-sframe
% tmp1 = fread(fp, [he_w he_h*mul], form, ofds(cnt), info.ByteOrder)';
%fileID , %dimensions of frame, %BitDepth, %offset, %ByteOrder
% imData{cnt}=cast(tmp1,form);
%end
%ieee-le.l64
imData=cell2mat(imData);
imData=reshape(imData,[he_h*mul,he_w,framenum]);%Nadav: Why is reshape needed if the values at each index of the cell array are already organized into hight by weight matrix.
fclose(fp);
display('Finished reading images')
elseif strcmpi(ext,'.hdf5') || strcmpi(ext,'.h5');
info = hdf5info(path_to_file);
dims = info.GroupHierarchy.Datasets.Dims;
if nargin < 2
sframe = 1;
end
if nargin < 3
num2read = dims(end)-sframe+1;
end
num2read = min(num2read,dims(end)-sframe+1);
imData = h5read(path_to_file,'/mov',[ones(1,length(dims)-1),sframe],[dims(1:end-1),num2read]);
else
error('Unknown file extension. Only .tiff and .hdf5 files are currently supported');
end