1function read(obj,filename,varargin)
2%READ Load sequence from file.
3% READ(seqObj, filename, ...) Read the given filename and load sequence
4% data into sequence object.
5%
6% optional parwameter 'detectRFuse' can be given to let the function
7% infer the currently missing flags concerning the intended use of the RF
8% pulses (excitation, refocusing, etc). These are important for the
9% k-space trajectory calculation
10%
11% Examples:
12% Load the sequence defined in gre.seq in my_sequences directory
13%
14% read(seqObj,'my_sequences/gre.seq')
15%
16% See also write
18detectRFuse=false;
19if ~isempty(varargin) && ~isempty(strfind(varargin{:},'detectRFuse'))
20 detectRFuse=true;
21end
23fid = fopen(filename);
25if fid<0
26 error('filed to open file ''%s''', filename);
27end
29% Clear sequence data
30%obj.blockEvents = [];
31obj.blockEvents = {};
32old_definitions = obj.definitions;
33obj.definitions = containers.Map();
34obj.gradLibrary = mr.EventLibrary();
35obj.shapeLibrary = mr.EventLibrary();
36obj.rfLibrary = mr.EventLibrary();
37obj.adcLibrary = mr.EventLibrary();
38%obj.delayLibrary = mr.EventLibrary();
39obj.trigLibrary = mr.EventLibrary();
40obj.labelsetLibrary = mr.EventLibrary();
41obj.labelincLibrary = mr.EventLibrary();
42obj.extensionStringIDs={};
43obj.extensionNumericIDs=[];
45version_combined=0;
46requiredDefs=struct('GradientRasterTime',false,'RadiofrequencyRasterTime',false,'AdcRasterTime',false,'BlockDurationRaster',false);
48% Load data from file
49while true
50 section = skipComments(fid);
51 if section == -1
52 break
53 end
55 switch section
56 case '[DEFINITIONS]'
57 obj.definitions = readDefinitions(fid);
58 v=obj.getDefinition('GradientRasterTime');
59 if ~isempty(v)
60 obj.gradRasterTime=v;
61 requiredDefs.GradientRasterTime=true;
62 end
63 v=obj.getDefinition('RadiofrequencyRasterTime');
64 if ~isempty(v)
65 obj.rfRasterTime=v;
66 requiredDefs.RadiofrequencyRasterTime=true;
67 end
68 v=obj.getDefinition('AdcRasterTime');
69 if ~isempty(v)
70 obj.adcRasterTime=v;
71 requiredDefs.AdcRasterTime=true;
72 end
73 v=obj.getDefinition('BlockDurationRaster');
74 if ~isempty(v)
75 obj.blockDurationRaster=v;
76 requiredDefs.BlockDurationRaster=true;
77 end
78 % octave doen't have struct2array(requiredDefs)
79 requiredDefsC=struct2cell(requiredDefs);
80 if version_combined >= 1004000 && ~all([requiredDefsC{:}])
81 fn=fieldnames(requiredDefs);
82 fn=fn([requiredDefsC{:}]==0);
83 error(['Required definitions ' sprintf('%s ',fn{:}) 'are missing in the file']);
84 end
85 case '[SIGNATURE]'
86 tmpSignDefs = readDefinitions(fid);
87 if isKey(tmpSignDefs,'Type')
88 obj.signatureType=tmpSignDefs('Type');
89 end
90 if isKey(tmpSignDefs,'Hash')
91 obj.signatureValue=tmpSignDefs('Hash');
92 obj.signatureFile='Text'; % we are reading a text file, so much is known for sure
93 end
94 case '[VERSION]'
95 [version_major, ...
96 version_minor, ...
97 version_revision] = readVersion(fid);
98 assert(version_major == obj.version_major, ...
99 'Unsupported version_major %d', version_major)
100 %
101 version_combined=1000000*version_major+1000*version_minor+version_revision;
102 %
103 if version_combined < 1002000
104 error('Unsupported version %d.%d.%d, only file format revision 1.2.0 and above are supported', version_major, version_minor, version_revision);
105 end
106 if version_combined < 1003001
107 warning('Loading older Pulseq format file (version %d.%d.%d) some code may function not as expected', version_major, version_minor, version_revision);
108 end
109 if version_combined >= 1005000 && detectRFuse
110 warning('Option ''detectRFuse'' is not supported for file format version 1.5.0 and above');
111 detectRFuse=false;
112 end
113 case '[BLOCKS]'
114 if ~exist('version_major')
115 error('Pulseq file MUST include [VERSION] section prior to [BLOCKS] section');
116 end
117 [obj.blockEvents,obj.blockDurations,delayInd_tmp] = readBlocks(fid, obj.blockDurationRaster, version_combined);
118 case '[RF]'
119 if version_combined >= 1005000
120 obj.rfLibrary = readEvents(fid, [1 1 1 1 1e-6 1e-6 1 1 1 1 NaN]); % this is 1.5.x format
121 elseif version_combined >= 1004000
122 obj.rfLibrary = readEvents(fid, [1 1 1 1 1e-6 1 1]); % this is 1.4.x format
123 % we fix it below
124 else
125 obj.rfLibrary = readEvents(fid, [1 1 1 1e-6 1 1]); % this is 1.3.x and below
126 % we will have to scan through the library later after all the shapes have been loaded
127 end
128 case '[GRADIENTS]'
129 if version_combined >= 1005000
130 obj.gradLibrary = readEvents(fid, [1 1 1 1 1 1e-6], 'g' ,obj.gradLibrary); % this is 1.5.x format
131 elseif version_combined >= 1004000
132 obj.gradLibrary = readEvents(fid, [1 1 1 1e-6], 'g' ,obj.gradLibrary); % this is 1.4.x format
133 else
134 obj.gradLibrary = readEvents(fid, [1 1 1e-6], 'g' ,obj.gradLibrary); % this is 1.3.x and below
135 end
136 case '[TRAP]'
137 obj.gradLibrary = readEvents(fid, [1 1e-6 1e-6 1e-6 1e-6], 't', obj.gradLibrary);
138 case '[ADC]'
139 if version_combined >= 1005000
140 obj.adcLibrary = readEvents(fid, [1 1e-9 1e-6 1 1 1 1 1]); % this is 1.5.x format
141 else
142 obj.adcLibrary=readEvents(fid, [1 1e-9 1e-6 1 1]); % this is 1.4.x and older format
143 % for now we don't have the phase vector in the ADC library
144 %obj.adcLibrary.data = [obj.adcLibrary.data(:,1:3) 0 obj.adcLibrary.data(:,4:5)]; % import from the older format
145 end
146 case '[DELAYS]'
147 if version_combined >= 1004000
148 error('Pulseq file revision 1.4.0 and above MUST NOT contain the [DELAYS] section');
149 end
150 tmp_delayLibrary = readEvents(fid, 1e-6);
151 case '[SHAPES]'
152 obj.shapeLibrary = readShapes(fid, (version_major==1 && version_minor<4));
153 case '[EXTENSIONS]'
154 obj.extensionLibrary = readEvents(fid);
155 otherwise
156 if strncmp('extension', section, 9)
157 extension=section(11:end);
158 if strncmp('TRIGGERS', extension, 8)
159 id=str2num(extension(9:end));
160 obj.setExtensionStringAndID('TRIGGERS',id);
161 obj.trigLibrary = readEvents(fid, [1 1 1e-6 1e-6]);
162 elseif strncmp('LABELSET', extension, 8)
163 id=str2num(extension(9:end));
164 obj.setExtensionStringAndID('LABELSET',id);
165 obj.labelsetLibrary = readAndParseEvents(fid,[],@str2num,@(s)find(strcmp(mr.getSupportedLabels,s)));
166 elseif strncmp('LABELINC', extension, 8)
167 id=str2num(extension(9:end));
168 obj.setExtensionStringAndID('LABELINC',id);
169 obj.labelincLibrary = readAndParseEvents(fid,[],@str2num,@(s)find(strcmp(mr.getSupportedLabels,s)));
170 elseif strncmp('DELAYS', extension, 6)
171 id=str2num(extension(7:end));
172 obj.setExtensionStringAndID('DELAYS',id);
173 obj.softDelayLibrary = readAndParseEvents(fid,[],@str2num,@(s) 1e-6*str2num(s),@str2num,@(s) parseSoftDelayHint(s, obj));
174 elseif strncmp('RF_SHIMS', extension, 8)
175 id=str2num(extension(9:end));
176 obj.setExtensionStringAndID('RF_SHIMS',id);
177 obj.rfShimLibrary = readAndParseEvents(fid,@preprocRfShimData);
178 elseif strncmp('ROTATIONS', extension, 9)
179 id=str2num(extension(10:end));
180 obj.setExtensionStringAndID('ROTATIONS',id);
181 obj.rotationLibrary = readEvents(fid);
182 for i=1:length(obj.rotationLibrary.data)
183 obj.rotationLibrary.data(i).array=mr.aux.quat.normalize(obj.rotationLibrary.data(i).array);
184 end
185 else
186 warning('Ignoring unknown extension, input string: %s', extension);
187 exts=regexp(extension, '(\s+)','split');
188 obj.setExtensionStringAndID(exts{1}, str2num(exts{2}));
189 skipSection(fid);
190 end
191 else
192 error('Unknown section code: %s', section);
193 end
194 end
195end
196fclose(fid);
198% fix sequence data imported from older verisons
199if version_combined < 1002000
200 error('Unsupported version %07d, only file format revision 1.2.0 (1002000) and above are supported', version_combined);
201end
203% a special case for ADCs as the format for them has only been updated once (in v1.5.0)
204% we have to do it first because seq.getBlock is used in the next version porting code section (version_combined < 1004000)
205if version_combined < 1005000
206 % scan though the ADCs and add empty phase shape IDs
207 for i=1:length(obj.adcLibrary.data)
208 obj.adcLibrary.update_data(...
209 obj.adcLibrary.keys(i), ...
210 obj.adcLibrary.data(i).array, ...
211 [obj.adcLibrary.data(i).array(1:3) 0 0 obj.adcLibrary.data(i).array(4:5) 0]); % add empty freqPPM, phasePPM and phase_id fields
212 end
213end
215% fix blocks, gradients and RF objects imported from older versions (< v1.4.0)
216if version_combined < 1004000
217 % fix definitions which are be missing in older files
218 if ~obj.definitions.isKey('GradientRasterTime')
219 obj.setDefinition('GradientRasterTime', obj.gradRasterTime);
220 end
221 if ~obj.definitions.isKey('RadiofrequencyRasterTime')
222 obj.setDefinition('RadiofrequencyRasterTime', obj.rfRasterTime);
223 end
224 if ~obj.definitions.isKey('AdcRasterTime')
225 obj.setDefinition('AdcRasterTime', obj.adcRasterTime);
226 end
227 if ~obj.definitions.isKey('BlockDurationRaster')
228 obj.setDefinition('BlockDurationRaster', obj.blockDurationRaster);
229 end
231 % scan through the RF objects
232 obj.rfLibrary.type(obj.rfLibrary.keys) = 'u'; % undefined for now, we'll attempt the type detection later (see below)
233 for i=1:length(obj.rfLibrary.data)
234 % % need to (partially) decode the magnitude shape to find out the pulse duration
235 %magSamples = obj.shapeLibrary.data(obj.rfLibrary.data(i).array(2)).array(1);
236 % % create time shape
237 %timeShape = mr.compressShape((1:magSamples)-0.5); % time shape is stored in units of RF raster
238 %data = [timeShape.num_samples timeShape.data];
239 %timeID = obj.shapeLibrary.find_or_insert(data);
240 rf=rmfield(obj.rfFromLibData([obj.rfLibrary.data(i).array(1:3) 0 0 obj.rfLibrary.data(i).array(4) 0 0 obj.rfLibrary.data(i).array(5:6)],'u'),'center');
241 center=mr.calcRfCenter(rf);
242 obj.rfLibrary.update_data(...
243 obj.rfLibrary.keys(i), ...
244 obj.rfLibrary.data(i).array, ...
245 [obj.rfLibrary.data(i).array(1:3) 0 center obj.rfLibrary.data(i).array(4) 0 0 obj.rfLibrary.data(i).array(5:6)]); % 0 between (4) and (5:6) are the freqPPM and phasePPM
246 end
248 % scan through the gradient objects and update 't'-s (trapezoids) und 'g'-s (free-shape gradients)
249 for i=1:length(obj.gradLibrary.data)
250 if obj.gradLibrary.type(i)=='t' % we need to fix some trapezoids, namely ones having zero amplitude and zero ramp times
251 if obj.gradLibrary.data(i).array(2)==0
252 if abs(obj.gradLibrary.data(i).array(1))==0 && obj.gradLibrary.data(i).array(3) > 0
253 obj.gradLibrary.update_data(...
254 obj.gradLibrary.keys(i), ...
255 obj.gradLibrary.data(i).array, ...
256 [obj.gradLibrary.data(i).array(1) obj.gradRasterTime obj.gradLibrary.data(i).array(3)-obj.gradRasterTime obj.gradLibrary.data(i).array(4:5)],...
257 obj.gradLibrary.type(i));
258 end
259 end
260 if obj.gradLibrary.data(i).array(4)==0
261 if abs(obj.gradLibrary.data(i).array(1))==0 && obj.gradLibrary.data(i).array(3) > 0
262 obj.gradLibrary.update_data(...
263 obj.gradLibrary.keys(i), ...
264 obj.gradLibrary.data(i).array, ...
265 [obj.gradLibrary.data(i).array(1:2) obj.gradLibrary.data(i).array(3)-obj.gradRasterTime obj.gradRasterTime obj.gradLibrary.data(i).array(5)],...
266 obj.gradLibrary.type(i));
267 end
268 end
269 end
270 if obj.gradLibrary.type(i)=='g'
271 % % need to (partially) decode the shape to find out the duration
272 %nSamples = obj.shapeLibrary.data(obj.gradLibrary.data(i).array(2)).array(1);
273 % % create time shape
274 %timeShape = mr.compressShape((1:nSamples)-0.5); % time shape is stored in units of grad raster
275 %data = [timeShape.num_samples timeShape.data];
276 %timeID = obj.shapeLibrary.find_or_insert(data);
277 obj.gradLibrary.update_data(...
278 obj.gradLibrary.keys(i), ...
279 obj.gradLibrary.data(i).array, ...
280 [obj.gradLibrary.data(i).array(1) NaN NaN obj.gradLibrary.data(i).array(2) 0 obj.gradLibrary.data(i).array(3)], ... % we use NaNs to label the non-initialized first/last fields. These will be restored in the code below
281 'g');
282 end
283 end
285 % for versions prior to 1.4.0 blockDurations have not been initialized
286 obj.blockDurations=zeros(1,length(obj.blockEvents));
287 % scan trhough blocks and calculate durations
288 for iB = 1:length(obj.blockEvents)
289 b=obj.getBlock(iB);
290 if delayInd_tmp(iB) > 0
291 b.delay.type = 'delay';
292 b.delay.delay = tmp_delayLibrary.data(delayInd_tmp(iB)).array;
293 end
294 obj.blockDurations(iB)=mr.calcDuration(b);
295 end
296elseif version_combined < 1005000
297 % port from v1.4.x : RF, ADC and GRAD objects need to be updated
298 % this needs to be done on the level of the libraries, because getBlock will fail
300 % scan though the RFs and add center, freqPPM, phasePPM and use fields
301 obj.rfLibrary.type(obj.rfLibrary.keys) = 'u'; % undefined for now, we'll attemp the type detection later (see below)
302 for i=1:length(obj.rfLibrary.data)
303 % use goes into the type field, and this is done separately
304 rf=rmfield(obj.rfFromLibData([obj.rfLibrary.data(i).array(1:4) 0 obj.rfLibrary.data(i).array(5) 0 0 obj.rfLibrary.data(i).array(6:7)],'u'),'center');
305 center=mr.calcRfCenter(rf);
306 obj.rfLibrary.update_data(...
307 obj.rfLibrary.keys(i), ...
308 obj.rfLibrary.data(i).array, ...
309 [obj.rfLibrary.data(i).array(1:4) center obj.rfLibrary.data(i).array(5) 0 0 obj.rfLibrary.data(i).array(6:7)]); % 0 between (5) and (6:7) are the freqPPM and phasePPM
310 end
311 % scan through the gradient objects and update 'g'-s (free-shape gradients)
312 for i=1:length(obj.gradLibrary.data)
313 if obj.gradLibrary.type(i)=='g'
314 obj.gradLibrary.update_data(...
315 obj.gradLibrary.keys(i), ...
316 obj.gradLibrary.data(i).array, ...
317 [obj.gradLibrary.data(i).array(1) NaN NaN obj.gradLibrary.data(i).array(2:4)], ... % we use NaNs to label the non-initialized first/last fields. These will be restored in the code below
318 'g');
319 end
320 end
321end
324% another run through for all older versions
325if version_combined < 1005000
326 gradChannels={'gx','gy','gz'};
327 gradPrevLast=zeros(1,length(gradChannels));
328 for iB = 1:length(obj.blockEvents)
329 b=obj.getBlock(iB);
330 block_duration=obj.blockDurations(iB);
331 %obj.blockDurations(iB)=block_duration;
332 % we also need to keep track of the event IDs because some Pulseq files written by external software may contain repeated entries so searching by content will fail
333 eventIDs=obj.blockEvents{iB};
334 processedGradIDs=zeros(1,length(gradChannels));
335 % update the objects by filling in the fields not contained in the
336 % pulseq file
337 for j=1:length(gradChannels)
338 grad=b.(gradChannels{j});
339 if isempty(grad)
340 gradPrevLast(j)=0;
341 continue;
342 end
343 if strcmp(grad.type,'grad')
344 if grad.delay>0
345 gradPrevLast(j)=0;
346 end
347 if isfield(grad,'first') && isfinite(grad.first)
348 continue;
349 end
350 grad.first = gradPrevLast(j);
351 % is this an extended trapezoid?
352 if grad.time_id~=0
353 grad.last=grad.waveform(end);
354 grad_duration=grad.delay+grad.tt(end);
355 else
356 % restore samples on the edges of the gradient raster intervals
357 % for that we need the first sample
358 odd_step1=[grad.first 2*grad.waveform'];
359 odd_step2=odd_step1.*(mod(1:length(odd_step1),2)*2-1);
360 waveform_odd_rest=(cumsum(odd_step2).*(mod(1:length(odd_step2),2)*2-1))';
361 grad.last = waveform_odd_rest(end);
362 grad_duration=grad.delay+length(grad.waveform)*obj.gradRasterTime;
363 end
364 % bookkeeping
365 gradPrevLast(j) = grad.last;
366 if grad_duration+eps<block_duration
367 gradPrevLast(j)=0;
368 end
369 %b.(gradChannels{j})=grad;
370 % update library object
371 % this does not work s we don't know how the amplitude was defined
372 % amplitude = max(abs(grad.waveform));
373 % if amplitude>0
374 % [~,~,fnz]=find(grad.waveform,1); % find the first non-zero value and make it positive
375 % amplitude=amplitude*sign(fnz);
376 % end
377 % need to recover the amplidute from the library data directly...
378 id=eventIDs(j+2);
379 if j>1 && any(processedGradIDs(1:j)==id)
380 continue; % avoid repeated updates if the same gradient is applied on differen gradient axes
381 end
382 processedGradIDs(j)=id;
383 amplitude=obj.gradLibrary.data(id).array(1);
384 %
385 old_data = [amplitude NaN NaN grad.shape_id grad.time_id grad.delay];
386 new_data = [amplitude grad.first grad.last grad.shape_id grad.time_id grad.delay];
387 update_data(obj.gradLibrary, id, old_data, new_data,'g');
388 else
389 gradPrevLast(j)=0;
390 end
391 end
392 %% copy updated objects back into the event library
393 %obj.setBlock(iB,b);
394 end
397%for iB=1:size(obj.blockEvents,1)
398% % update the objects by filling in the fields not contained in the
399% % pulseq file
400% for j=1:length(gradChannels)
401% grad=b.(gradChannels{j});
402% if isempty(grad)
403% continue;
404% end
405% if strcmp(grad.type,'grad')
406% grad.first = grad.waveform(1); % MZ: eventually we should use extrapolation by 1/2 gradient rasters here
407% grad.last = grad.waveform(end);
408% b.(gradChannels{j})=grad;
409% end;
410% end
411% % copy updated objects back into the event library
412% obj.setBlock(iB,b);
413end
415if detectRFuse
416 % find the RF pulses, list flip angles
417 % and work around the current (rev 1.2.0) Pulseq file format limitation
418 % that the RF pulse use is not stored in the file
419 for k=obj.rfLibrary.keys
420 libData=obj.rfLibrary.data(k).array;
421 rf=obj.rfFromLibData(libData,'u');
422 %flipAngleDeg=abs(sum(rf.signal))*rf.t(1)*360; %we use rfex.t(1) in place of opt.system.rfRasterTime
423 flipAngleDeg=abs(sum(rf.signal(1:end-1).*(rf.t(2:end)-rf.t(1:end-1))))*360;
424 offresonance_ppm=1e6*rf.freqOffset/obj.sys.B0/obj.sys.gamma;
425 % fix library %%%% if length(obj.rfLibrary.type)>=eventInd(2)
426 if flipAngleDeg < 90.01 % we add 0.01 degree to account for rounding errors which we've experienced for very short RF pulses
427 obj.rfLibrary.type(k) = 'e';
428 else
429 if rf.shape_dur > 6e-3 && offresonance_ppm >= -3.5 && offresonance_ppm <= -3.4 % approx -3.45 ppm
430 obj.rfLibrary.type(k) = 's'; % saturation (fat-sat)
431 else
432 obj.rfLibrary.type(k) = 'r';
433 end
434 end
435% % fix libData
436% if length(libData) < 9
437% if flipAngleDeg < 90.01 % we add 0.01 degree to account for rounding errors which we've experienced for very short RF pulses
438% libData(9) = 0; % or 1 ?
439% else
440% libData(9) = 2; % or 1 ?
441% end
442% obj.rfLibrary.data(k).array=libData;
443% end
444 end
445end
447return
449%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
450%%%%%%%%%%%%%%%%%%%%%%% Helper functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
452 function def = readDefinitions(fid)
453 %readDefinitions Read the [DEFINITIONS] section of a sequence file.
454 % defs=readDefinitions(fid) Read user definitions from file
455 % identifier of an open MR sequence file and return a map of
456 % key/value entries.
458 def = containers.Map;
459 %line = localStrip(fgetl(fid));
460 line = localStrip(skipComments(fid));
461 while ischar(line) && ~(isempty(line) || line(1) == '#')
462 tok = textscan(line, '%s');
463 def(tok{1}{1}) = str2double(tok{1}(2:end));
464 if ~all(isfinite(def(tok{1}{1})))
465 def(tok{1}{1}) = strtrim(line((length(tok{1}{1})+2):end));
466 end
467 line = fgetl(fid);
468 end
469 end
471 function [major, minor, revision] = readVersion(fid)
472 %readVersion Read the [VERSION] section of a sequence file.
473 % defs=readVersion(fid) Read Pulseq version from file
474 % identifier of an open MR sequence file and return it
476 major = [];
477 minor = [];
478 revision = [];
479 line = fgetl(fid);
480 while ischar(line) && ~(isempty(line) || line(1)=='#')
481 tok = textscan(line,'%s');
482 switch tok{1}{1}
483 case 'major'
484 major = str2double(tok{1}(2:end));
485 case 'minor'
486 minor = str2double(tok{1}(2:end));
487 case 'revision'
488 revision = str2double(tok{1}(2:end));
489 end
490 line = fgetl(fid);
491 end
492 end
494 function [eventTable,blockDurations,delayIDs_tmp] = readBlocks(fid, blockDurationRaster, version_combined)
495 %readBlocks Read the [BLOCKS] section of a sequence file.
496 % library=readBlocks(fid) Read blocks from file identifier of an
497 % open MR sequence file and return the event table.
499 eventTable = {};
500 blockDurations = [];
501 delayIDs_tmp = [];
502 line = fgetl(fid);
503 while ischar(line) && ~(isempty(line) || line(1) == '#')
504 blockEvents = sscanf(line, '%f')';
505 %eventTable = [eventTable; blockEvents(2:end)];
506 if version_combined<=1002001
507 eventTable{blockEvents(1)} = [0 blockEvents(3:end) 0];
508 else
509 eventTable{blockEvents(1)} = [0 blockEvents(3:end)];
510 end
511 if version_combined>=1004000
512 blockDurations(blockEvents(1)) = blockEvents(2)*blockDurationRaster;
513 else
514 delayIDs_tmp(blockEvents(1)) = blockEvents(2);
515 end
516 line = fgetl(fid);
517 end
518 end
520 function str = format_helper(scale)
521 if isfinite(scale)
522 str='%f ';
523 else
524 str='%c ';
525 end
526 end
528 function eventLibrary = readEvents(fid, scale, type, eventLibrary)
529 %readEvents Read an event section of a sequence file.
530 % library=readEvents(fid) Read event data from file identifier of
531 % an open MR sequence file and return a library of events.
532 %
533 % library=readEvents(fid,scale) Read event data and scale
534 % elements according to column vector scale.
535 %
536 % library=readEvents(fid,scale,type) Attach the type string to
537 % elements of the library.
538 %
539 % library=readEvents(...,library) Append new events to the given
540 % library.
541 if nargin < 2
542 scale = 1;
543 format='%f';
544 type_idx=[];
545 data_mask=[];
546 else
547 % new in v1.5.0 : generate format string; NaN labels character param(s)
548 format=['%f ' cell2mat(arrayfun(@format_helper,scale,'UniformOutput',false))];
549 format(end)=[]; % matlab is so incredibly ugly!
550 data_mask=isfinite(scale);
551 type_idx=find(~data_mask);
552 if length(type_idx)>2
553 error('Only one type field (marked as NaN) can be provided');
554 end
555 if isempty(type_idx)
556 data_mask=[];
557 end
558 end
559 if nargin < 4
560 eventLibrary = mr.EventLibrary();
561 end
562 %
563 line = fgetl(fid);
564 while ischar(line) && ~(isempty(line) || line(1) == '#')
565 data = sscanf(line,format)';
566 id = data(1);
567 if ~isempty(type_idx)
568 type=char(data(type_idx+1)); % need +1 because of the eventID in the first position
569 end
570 data = scale.*data(2:end);
571 if nargin < 3 && isempty(type_idx)
572 if isempty(data_mask)
573 eventLibrary.insert(id, data);
574 else
575 eventLibrary.insert(id, data(data_mask));
576 end
577 else
578 if isempty(data_mask)
579 eventLibrary.insert(id, data, type);
580 else
581 eventLibrary.insert(id, data(data_mask), type);
582 end
583 end
585 line=fgetl(fid);
586 end
587 end
589 function eventLibrary = readAndParseEvents(fid, preproc, varargin)
590 %readAndParseEvents Read an event section of a sequence file.
591 % library=readAndParseEvents(fid) Read event data from file
592 % identifier of an open MR sequence file and return a library of
593 % events.
594 %
595 % library=readAndParseEvents(fid,[],parser1,parser2,...) Read event
596 % data and convert the elements using to the provided parser.
597 % Default parser is str2num()
598 %
599 % library=readAndParseEvents(fid,preproc,...) Read event data
600 % from file and apply the preproc() to the acquired data
601 % line-by-line prior to adding them to the event library
602 %
604 eventLibrary = mr.EventLibrary();
605 line = fgetl(fid);
606 while ischar(line) && ~(isempty(line) || line(1) == '#')
607 datas=regexp(line, '(\s+)','split');
608 data=zeros(1,length(datas)-1);
609 id = str2num(datas{1});
610 for i=2:length(datas)
611 if i>nargin-1
612 data(i-1) = str2num(datas{i});
613 else
614 data(i-1) = varargin{i-1}(datas{i});
615 end
616 end
618 if ~exist('preproc','var') || isempty(preproc)
619 eventLibrary.insert(id, data);
620 else
621 eventLibrary.insert(id, preproc(data));
622 end
624 line=fgetl(fid);
625 end
626 end
628 function skipSection(fid)
629 %skipSection Read an event section of a sequence file without
630 % interpreting it.
631 %
632 line = fgetl(fid);
633 while ischar(line) && ~(isempty(line) || line(1) == '#')
634 line=fgetl(fid);
635 end
636 end
638 function shapeLibrary = readShapes(fid, forceConvertUncompressed)
639 %readShapes Read the [SHAPES] section of a sequence file.
640 % library=readShapes(fid) Read shapes from file identifier of an
641 % open MR sequence file and return a library of shapes.
643 shapeLibrary = mr.EventLibrary();
644 line = skipComments(fid);
645 while ~(~ischar(line) || isempty(line) || ~strcmp(line(1:8), 'shape_id'))
646 tok = textscan(line, '%s');
647 id = str2double(tok{1}(2));
648 line = skipComments(fid);
649 tok = textscan(line, '%s');
650 num_samples = str2double(tok{1}(2));
651 data = [];
652 line = skipComments(fid); % first sample
653 while ischar(line) && ~(isempty(line) || line(1) == '#')
654 data = [data sscanf(line, '%f')];
655 %data = [data single(sscanf(line, '%f'))]; % C-code uses single precision and we had problems already due to the rounding during reading in of the shapes...
656 line = fgetl(fid);
657 end
659 line = skipComments(fid, true); % MZ: second parameter to prevent readShapes from reading into the next section (long-standing bug)
661 % check if conversion is needed: in v1.4.x we use length(data)==num_samples
662 % as a marker for the uncompressed (stored) data. In older versions this condition could occur by chance
663 if forceConvertUncompressed && length(data)==num_samples
664 shape.data=data;
665 shape.num_samples=num_samples;
666 shape = mr.compressShape(mr.decompressShape(shape,true));
667 data = [shape.num_samples shape.data];
668 else
669 data = [num_samples data];
670 end
671 shapeLibrary.insert(id, data);
672 end
673 end
675 function nextLine = skipComments(fid, stopBeforeSection)
676 %skipComments Read lines of skipping blank lines and comments.
677 % line=skipComments(fid) Read lines from valid file identifer and
678 % return the next non-comment line.
680 if (nargin<2)
681 stopBeforeSection=false;
682 end
684 tmpPos=ftell(fid);
685 line = fgetl(fid);
686 while ischar(line) && (isempty(line) || line(1) == '#')
687 tmpPos=ftell(fid);
688 line = fgetl(fid);
689 end
690 if ischar(line)
691 if stopBeforeSection && line(1)=='['
692 fseek(fid,tmpPos,-1); % restore the file position
693 nextLine = ''; % feasible (non-error) dummy return
694 else
695 nextLine = line;
696 end
697 else
698 nextLine = -1;
699 end
700 end
702 function id=parseSoftDelayHint(s, seq)
703 try
704 id=seq.softDelayHints1(s);
705 catch
706 id=seq.softDelayHints1.length()+1;
707 seq.softDelayHints1(s)=id;
708 seq.softDelayHints2{id}=s;
709 end
710 end
712 function data_out=preprocRfShimData(data)
713 if length(data)~=data(1)*2+1
714 error('Error reading RF shim extension data');
715 end
716 data_out=data(2:end);
717 end
719 % for compatibility with Octave which has no strip()
720 function s=localStrip(s)
721 a=1;
722 b=length(s);
723 while a<=b && isspace(s(a)), a=a+1; end
724 while a<=b && isspace(s(b)), b=b-1; end
725 s=s(a:b);
726 end
727end