1classdef Sequence < handle
2 % Sequence Generate sequences and read/write sequence files.
3 % This class defines properties and methods to define a complete
4 % MR sequence including RF pulses, gradients, ADC events, etc.
5 %
6 % The class provides an implementation of the open MR sequence format
7 % defined by the Pulseq project.
8 % See http://pulseq.github.io/
9 %
10 % Sequence Properties:
11 % definitions - A list of custom definitions
12 %
13 % Sequence Methods:
14 % read - Load sequence from open MR sequence format
15 % write - Write sequence to open MR sequence format
16 %
17 % Sequence Static Methods:
18 % makeTrapezoid - Create a trapezoid gradient structure
19 %
20 % Examples:
21 %
22 % To read a sequence from file:
23 % read(seqObj,'my_sequences/gre.seq');
24 %
25 % To plot a sequence:
26 % plot(seqObj)
27 %
28 % See also demoRead.m, demoWrite.m
29 % Examples defining an MRI sequence and reading/writing files
30 %
31 % Kelvin Layton <kelvin.layton@uniklinik-freiburg.de>
32 % Maxim Zaitsev <maxim.zaitsev@uniklinik-freiburg.de>
34 % Private properties
35 %
36 properties(GetAccess = public, SetAccess = private)
37 version_major;
38 version_minor;
39 version_revision;
40 rfRasterTime; % RF raster time (system dependent)
41 gradRasterTime; % Gradient raster time (system dependent)
42 adcRasterTime; % minimum unit/increment of the ADC dwell time (system dependent)
43 blockDurationRaster; % unit/increment of the block duration (system dependent)
44 definitions % Optional sequence definitions
46 blockEvents; % Event table (references to events)
47 blockDurations; % List of block durations
48 rfLibrary; % Library of RF events
49 gradLibrary; % Library of gradient events
50 adcLibrary; % Library of ADC readouts
51 trigLibrary; % Library of trigger events ( referenced from the extentions library )
52 labelsetLibrary; % Library of Label(set) events ( reference from the extensions library )
53 labelincLibrary; % Library of Label(inc) events ( reference from the extensions library )
54 extensionLibrary; % Library of extension events. Extension events form single-linked zero-terminated lists
55 shapeLibrary; % Library of compressed shapes
56 rfShimLibrary; % Library of RF shimming events
57 softDelayLibrary; % Library of 'soft delay' extension events.
58 softDelayHints1; % Map of string hints that are the part of the 'soft delay' extension objects
59 softDelayHints2; % cell array of string hints that are the part of the 'soft delay' extension objects
60 rotationLibrary; % Library of the rotation extension objects (rotation quaternions)
61 extensionStringIDs; % string IDs of the used extensions (cell array)
62 extensionNumericIDs; % numeric IDs of the used extensions (numeric array)
64 gradCheckData; % struct caching date used for checking of extended gradients cthat cross block boundaries
66 signatureType; % type of the hashing function used, currently 'md5'
67 signatureFile; % which data were hashed, currently 'text' or 'bin' (used file format of the save function)
68 signatureValue; % the hash of the exported Pulse sequence
70 rfID2NameMap; % optional names of objects in the plot
71 adcID2NameMap; % optional names of objects in the plot
72 gradID2NameMap; % optional names of objects in the plot
74 % TRID helper (name -> numeric ID mapping)
75 tridName2Id; % containers.Map('char' -> int32)
76 tridId2Name; % cell array, index = numeric ID
77 tridHistory; % cell array storing TRID calls in order
79 sys;
80 end
82 methods
84 function obj = Sequence(varargin)
85 [obj.version_major, obj.version_minor, obj.version_revision] = mr.aux.version();
86 % version minor 3 will now support control events (8th column in the event table) mv4 supports/expects timing vectors for arbitrary grads
87 obj.definitions = containers.Map();
88 obj.gradLibrary = mr.EventLibrary();
89 obj.shapeLibrary = mr.EventLibrary();
90 obj.rfLibrary = mr.EventLibrary();
91 obj.adcLibrary = mr.EventLibrary();
92 obj.trigLibrary = mr.EventLibrary();
93 obj.labelsetLibrary = mr.EventLibrary();
94 obj.labelincLibrary = mr.EventLibrary();
95 obj.rfShimLibrary = mr.EventLibrary();
96 obj.rotationLibrary = mr.EventLibrary();
97 obj.extensionLibrary = mr.EventLibrary();
98 obj.extensionStringIDs={};
99 obj.extensionNumericIDs=[];
100 obj.softDelayLibrary = mr.EventLibrary();
101 obj.softDelayHints1 = containers.Map();
102 obj.softDelayHints2 = {};
103 obj.blockEvents = {};
105 if nargin<1
106 sys=mr.opts();
107 else
108 sys=varargin{1};
109 end
110 if ~isfield(sys, 'flag_trid') || isempty(sys.flag_trid)
111 sys.flag_trid = true;
112 end
113 sys.flag_trid = logical(sys.flag_trid);
115 obj.sys = sys;
116 obj.rfRasterTime = sys.rfRasterTime;
117 obj.gradRasterTime = sys.gradRasterTime;
118 obj.adcRasterTime = sys.adcRasterTime;
119 obj.blockDurationRaster = sys.blockDurationRaster;
120 obj.setDefinition('GradientRasterTime', obj.gradRasterTime);
121 obj.setDefinition('RadiofrequencyRasterTime', obj.rfRasterTime);
122 obj.setDefinition('AdcRasterTime', obj.adcRasterTime);
123 obj.setDefinition('BlockDurationRaster', obj.blockDurationRaster);
124 obj.signatureType='';
125 obj.signatureFile='';
126 obj.signatureValue='';
127 obj.rfID2NameMap = containers.Map('KeyType', 'int32', 'ValueType', 'char');
128 obj.adcID2NameMap = containers.Map('KeyType', 'int32', 'ValueType', 'char');
129 obj.gradID2NameMap = containers.Map('KeyType', 'int32', 'ValueType', 'char');
130 obj.tridName2Id = containers.Map('KeyType','char','ValueType','int32');
131 obj.tridId2Name = {};
132 obj.tridHistory = {};
133 obj.gradCheckData=struct('validForBlockNum',0,'lastGradVals', [0 0 0]);
135 end
137 function copyDefinitions(obj, otherSeq)
138 % copy all definitions from another sequence
139 % in future we may add optional include or exclude filters as
140 % optional parameters
141 obj.definitions=otherSeq.definitions;
142 end
145 function addTRID(obj, label_name)
146 %addTRID Add a GE TRID segment label (by name).
147 % The TRID label is ignored if obj.sys.flag_trid==false.
148 % label_name is mapped to a numeric TRID ID automatically
149 % (first occurrence defines the ID).
150 if ~isfield(obj.sys,'flag_trid') || ~obj.sys.flag_trid
151 return;
152 end
153 id = obj.getOrCreateTridId(label_name);
154 obj.addBlock(mr.makeLabel('SET','TRID', double(id)));
155 end
157 % See read.m
158 read(obj,filename,varargin)
160 % See write.m
161 write(obj,filename,create_signature)
163 % See write_v141.m
164 write_v141(obj,filename,create_signature)
166 % See write.m
167 write_file(obj,filename)
169 % See readBinary.m
170 readBinary(obj,filename);
172 % See writeBinary.m
173 writeBinary(obj,filename);
176 % See calcPNS.m
177 [ok, pns_norm, pns_comp, t_axis]=calcPNS(obj,hardware,doPlots,calcCNS)
179 %see calcMomentsBtensor.m
180 [B, m1, m2, m3] = calcMomentsBtensor(obj, calcB, calcm1, calcm2, Ndummy, calcm3)
182 % See testReport.m
183 [ report ] = testReport( obj, varargin )
185 % See autoLabel.m
186 [labels, aux] = autoLabel(seq, varargin)
188 % See gradSpectrum.m
189 [R, Rax, F] = gradSpectrum(obj, FB, fmax, plt)
191 function [duration, numBlocks, eventCount]=duration(obj)
192 % duration()
193 % Returns the total duration of the sequence
194 % optionally returns the total count of events
195 %
197 % Loop over blocks and gather statistics
198 numBlocks = length(obj.blockEvents);
199 if numBlocks>0 && nargout>2
200 eventCount=zeros(size(obj.blockEvents{1}));
201 end
202 duration=0;
203 for iB=1:numBlocks
204 if nargout>2
205 eventCount = eventCount + (obj.blockEvents{iB}>0);
206 end
207 duration=duration+obj.blockDurations(iB);
208 end
209 end
211 function [is_ok, errorReport]=checkTiming(obj)
212 % checkTiming()
213 % Checks timing (and some other parameters) of all blocks
214 % and objects in the sequence optionally returns a detailed
215 % error log as cell array of strings. This function also
216 % modifies the sequence object by adding the field
217 % "TotalDuration" to sequence definitions
218 %
220 % Loop over blocks and gather statistics
221 numBlocks = length(obj.blockEvents);
222 is_ok=true;
223 errorReport={};
224 totalDuration=0;
225 gradBook=struct();
226 for iB=1:numBlocks
227 b=obj.getBlock(iB);
228 % assemble cell array of events
229 %ev={b.rf, b.gx, b.gy, b.gz, b.adc, b.delay, b.ext};
230 %ind=~cellfun(@isempty,ev);
231 % the above does not work for ext because it may be
232 % missing from some blocks and may have multiple entries in
233 % others.
234 ind=~structfun(@isempty,b);
235 fn=fieldnames(b);
236 ev=cellfun(@(f) b.(f), fn(ind), 'UniformOutput', false);
237 [res, rep, dur] = mr.checkTiming(obj.sys,ev{:}); %ev{ind});
239 is_ok = (is_ok && res);
241 % check the stored block duration
242 if abs(dur-obj.blockDurations(iB))>eps
243 rep = [rep ' inconsistency between the stored block duration and the duration of the block content'];
244 is_ok = false;
245 dur=obj.blockDurations(iB);
246 end
248 % check that block duration is aligned to the blockDurationRaster
249 bd=obj.blockDurations(iB)/obj.blockDurationRaster;
250 bdr=round(bd);
251 if abs(bdr-bd)>=1e-6
252 rep = [rep ' block duration is not aligned to the blockDurationRaster'];
253 is_ok = false;
254 end
256 % check RF dead times
257 if ~isempty(b.rf)
258 if b.rf.delay-b.rf.deadTime < -eps
259 rep = [rep ' delay of ' num2str(b.rf.delay*1e6) 'us is smaller than the RF dead time ' num2str(b.rf.deadTime*1e6) 'us'];
260 is_ok = false;
261 end
262 if b.rf.delay+b.rf.t(end)+b.rf.ringdownTime-dur > eps
263 rep = [rep ' time between the end of the RF pulse at ' num2str((b.rf.delay+b.rf.t(end))*1e6) ' and the end of the block at ' num2str(dur*1e6) 'us is shorter than rfRingdownTime'];
264 is_ok = false;
265 end
266 if abs(b.rf.freqOffset) > obj.sys.maxFreqOffset || abs(b.rf.freqPPM*1e-6*obj.sys.gamma) > obj.sys.maxFreqOffset || abs(b.rf.freqOffset + b.rf.freqPPM*1e-6*obj.sys.gamma) > obj.sys.maxFreqOffset
267 rep = [rep ' frequency offset of the RF pulse exceeds the maximum allowed value of ' num2str(obj.sys.maxFreqOffset) 'Hz'];
268 is_ok = false;
269 end
270 end
272 % check ADC dead times, dwell times and numbers of samples
273 if ~isempty(b.adc)
274 if b.adc.delay-obj.sys.adcDeadTime < -eps
275 rep = [rep ' adc.delay<system.adcDeadTime'];
276 is_ok=false;
277 end
278 if b.adc.delay+b.adc.numSamples*b.adc.dwell+obj.sys.adcDeadTime-dur > eps
279 rep = [rep ' adc: system.adcDeadTime (post-adc) violation'];
280 is_ok=false;
281 end
282 if abs(b.adc.dwell/obj.sys.adcRasterTime-round(b.adc.dwell/obj.sys.adcRasterTime)) > 1e-10 % the check against eps was too strict
283 rep = [rep ' adc: dwell time is not an integer multiple of sys.adcRasterTime'];
284 is_ok=false;
285 end
286 if abs(b.adc.numSamples/obj.sys.adcSamplesDivisor-round(b.adc.numSamples/obj.sys.adcSamplesDivisor)) > eps
287 rep = [rep ' adc: numSamples is not an integer multiple of sys.adcSamplesDivisor'];
288 is_ok=false;
289 end
290 if abs(b.adc.freqOffset) > obj.sys.maxFreqOffset || abs(b.adc.freqPPM*1e-6*obj.sys.gamma) > obj.sys.maxFreqOffset || abs(b.adc.freqOffset + b.adc.freqPPM*1e-6*obj.sys.gamma) > obj.sys.maxFreqOffset
291 rep = [rep ' frequency offset of the ADC object exceeds the maximum allowed value of ' num2str(obj.sys.maxFreqOffset) 'Hz'];
292 is_ok = false;
293 end
294 end
296 % update report
297 if ~isempty(rep)
298 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' ' rep '\n' ] };
299 end
301 % check shaped gradients that may potentially end/start at non-zero values
302 gradBookCurr=struct();
303 if ~isempty(ev) && iscell(ev)
304 for en=1:length(ev)
305 if length(ev{en})==1 && isstruct(ev{en}) && strcmp(ev{en}.type,'grad') % length(ev{en})==1 excludes arrays of extensions
306 g=ev{en};
307 if g.first~=0
308 if g.delay~=0
309 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' ' g.channel ' gradient starts at a non-zero value but defines a delay\n' ] };
310 is_ok=false;
311 end
312 if ~isfield(gradBook, g.channel) || abs(gradBook.(g.channel)-g.first) > 1 % 1 Hz/m is ~23 nT/m - this tolerance originates from the library compression; otherwise we have to increase the number of digits when generating search strings in the library code... %1e-6 % todo: real physical tolarance for gradient amplitudes
313 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' ' g.channel ' gradient''s start value ' num2str(g.first) ' differs from the previous block end value\n' ] };
314 is_ok=false;
315 else
316 gradBook.(g.channel)=0; % reset as properly consumed
317 end
318 end
319 if abs(g.last)>eps
320 if abs(g.delay+g.shape_dur - dur) > eps
321 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' ' g.channel ' gradient ends at a non-zero value but does not last until the end of the block\n' ] };
322 is_ok=false;
323 end
324 gradBookCurr.(g.channel)=g.last; % update bookkeeping
325 end
326 end
327 end
328 end
330 % check soft delays
331 if isfield(b, 'softDelay') && ~isempty(b.softDelay)
332 if ~exist('softDelayState','var')
333 softDelayState={};
334 end
335 if b.softDelay.factor==0
336 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' soft delay ' b.softDelay.hint '/' num2str(b.softDelay.num) ' has factor parameter of 0 which is invalid\n' ] };
337 is_ok=false;
338 end
339 % calculate the default delay value based on the current block duration
340 def_del=(obj.blockDurations(iB)-b.softDelay.offset)*b.softDelay.factor;
341 if (b.softDelay.num>=0)
342 % remember or check for consistency
343 if length(softDelayState)<b.softDelay.num+1 || isempty(softDelayState{b.softDelay.num+1})
344 softDelayState{b.softDelay.num+1}=struct('def',def_del,'hint',b.softDelay.hint, 'blk', iB);
345 else
346 if abs(def_del-softDelayState{b.softDelay.num+1}.def)>1e-7 % what is the reasonable threshold?
347 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' soft delay ' b.softDelay.hint '/' num2str(b.softDelay.num) ': default duration derived from this block (' num2str(def_del*1e6) 'us) is inconsistent with the previous default (' num2str(softDelayState{b.softDelay.num+1}.def*1e6) 'us) that was derived from block ' num2str(softDelayState{b.softDelay.num+1}.blk) '\n' ] };
348 is_ok=false;
349 end
350 if ~strcmp(b.softDelay.hint, softDelayState{b.softDelay.num+1}.hint)
351 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' soft delay ' b.softDelay.hint '/' num2str(b.softDelay.num) ': soft delays with the same numeric ID are expected to share the same text hint but previous hint recorded in block ' num2str(softDelayState{b.softDelay.num+1}.blk) ' is ' softDelayState{b.softDelay.num+1}.hint '\n' ] };
352 is_ok=false;
353 end
354 end
355 else
356 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' contains a soft delay ' b.softDelay.hint ' with an invalid numeric ID' num2str(b.softDelay.num) '\n' ] };
357 is_ok=false;
358 end
359 end
361 % check whether all gradient bookkeeping values have been properly consumed
362 if dur~=0
363 % octave has no struct2array()
364 gradBookC = struct2cell(gradBook);
365 if any(0~=[gradBookC{:}])
366 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' some gradients in the previous non-empty block are ending at non-zero values but are not continued here\n' ] };
367 is_ok=false;
368 end
369 gradBook=gradBookCurr;
370 end
372 %
373 totalDuration = totalDuration+dur;
374 end
376 % check whether all gradients in the last block are ramped down properly
377 if ~isempty(ev) && iscell(ev)
378 for en=1:length(ev)
379 if length(ev{en})==1 && isstruct(ev{en}) && strcmp(ev{en}.type,'grad') % length(ev{en})==1 excludes arrays of extensions
380 if ev{en}.last~=0 % must be > sys.slewRate*sys.gradRasterTime
381 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' gradients do not ramp to 0 at the end of the sequence\n' ] };
382 is_ok=false;
383 end
384 end
385 end
386 end
388 prevTotalDuration=obj.getDefinition('TotalDuration');
389 if ~isempty(prevTotalDuration) && abs(prevTotalDuration-totalDuration)>1e-9
390 errorReport = { errorReport{:}, [ ' TotalDuration definition of ' sprintf('%.9g', prevTotalDuration) 's was present in the sequence, but was incorrect. It is now ' sprintf('%.9g', totalDuration) 's\n' ] };
391 is_ok=false;
392 end
393 obj.setDefinition('TotalDuration', totalDuration);
394 end
396 function value=getDefinition(obj,key)
397 %getDefinition Return the values of custom definition.
398 % val=getDefinitions(seqObj,key) Return value of the
399 % definition specified by the key.
401 % These definitions can be added manually or read from the
402 % header of a sequence file defined in the sequence header.
403 % An empty array is return if the key is not defined.
404 %
405 % See also setDefinition
406 if isKey(obj.definitions,key)
407 value = obj.definitions(key);
408 else
409 value = [];
410 end
411 end
413 function setDefinition(seqObj,key,val)
414 %setDefinition Modify a custom definition of the sequence.
415 % setDefinition(seqObj,def,val) Set the user definition 'key'
416 % to value 'val'. If the definition does not exist it will be
417 % created.
418 %
419 % See also getDefinition
420 if strcmp(key,'FOV')
421 % issue a warning if FOV is too large e.g. is in mm
422 if max(val)>1
423 warning('WARNING: definition FOV uses values exceeding 1m. New Pulseq interpreters expect values in units of meters!\n');
424 end
425 end
426 seqObj.definitions(key)=val;
427 end
429 function addBlock(obj,varargin)
430 %addBlock Add a new block to the sequence.
431 % addBlock(obj, blockStruct) Adds a sequence block with
432 % provided as a block struture
433 %
434 % addBlock(obj, e1, e2, ...) Adds a block with multiple
435 % events e1, e2, etc.
436 %
437 % addBlock(obj, duration, e1, e2, ...) Create a new block
438 % with the given predefined duration populated with events
439 % e1, e2, etc. If the duration of any of the events exceeds
440 % the desired duration an error will be thrown.
441 %
442 % See also setBlock, makeAdc, makeTrapezoid, makeSincPulse
443 %setBlock(obj,size(obj.blockEvents,1)+1,varargin{:});
444 setBlock(obj,length(obj.blockEvents)+1,varargin{:});
445 end
447 function iB=findBlockByTime(obj,t)
448 if nargin<3
449 nonEmpty=true;
450 end
451 iB=find(cumsum(obj.blockDurations)>t,1);
452 if iB>length(obj.blockDurations);
453 iB=[]; % or length(obj.blockDurations)
454 end
455 assert(obj.blockDurations(iB)>0);
456 %if nonEmpty && ~isempty(iB)
457 % iB=find(obj.blockDurations(1:iB)>0,1,'last');
458 %end
459 end
461 function modGradAxis(obj,axis,modifier)
462 %modGradAxis Invert or scale all gradinents along the corresponding
463 % axis/channel. The function acts on all gradient objects
464 % already added to the sequence object
465 %
466 channelNum = find(strcmp(axis, ...
467 {'x', 'y', 'z'}));
468 otherChans = find(~strcmp(axis, ...
469 {'x', 'y', 'z'}));
470 % go through all event table entries and list gradient
471 % objects in the library
472 %paren = @(x, varargin) x(varargin{:}); % anonymous function to access the array on the fly
473 paren2 = @(x, varargin) x(:,varargin{:}); % anonymous function to access the array on the fly
474 %allGradEvents = paren(vertcat(obj.blockEvents{:}),:,3:5);
475 allGradEvents = paren2(vertcat(obj.blockEvents{:}),3:5);
477 selectedEvents=unique(allGradEvents(:,channelNum));
478 selectedEvents=selectedEvents(0~=selectedEvents); % elliminate 0
479 otherEvents=unique(allGradEvents(:,otherChans));
480 assert(isempty(intersect(selectedEvents,otherEvents)),'ERROR: the same gradient event is used on multiple axes, this is not yet supported by modGradAxis()');
482 for i = 1:length(selectedEvents)
483 %type = obj.gradLibrary.type(i);
484 %libData = obj.gradLibrary.data(i).array;
485 %if strcmp(grad.type,'grad')
486 % amplitude = libData(1);
487 %else
488 % %grad.amplitude = libData(1);
489 %end
490 %
491 % based on the above we just patch the first element of the
492 % gradient library data entries
493 obj.gradLibrary.data(selectedEvents(i)).array(1)=modifier*obj.gradLibrary.data(selectedEvents(i)).array(1);
494 if obj.gradLibrary.type(selectedEvents(i))=='g' && obj.gradLibrary.lengths(selectedEvents(i))==5
495 % need to update .first and .last fields
496 obj.gradLibrary.data(selectedEvents(i)).array(2)=modifier*obj.gradLibrary.data(selectedEvents(i)).array(2); % change in v150
497 obj.gradLibrary.data(selectedEvents(i)).array(3)=modifier*obj.gradLibrary.data(selectedEvents(i)).array(3); % change in v150
498 end
499 end
500 end
502 function flipGradAxis(obj, axis)
503 %flipGradAxis Invert all gradinents along the corresponding
504 % axis/channel. The function acts on all gradient objects
505 % already added to the sequence object
506 %
507 modGradAxis(obj,axis,-1);
508 end
510 function rf = rfFromLibData(obj, libData, use)
511 rf.type = 'rf';
513 amplitude = libData(1);
514 magShape = libData(2);
515 phaseShape = libData(3);
516 shapeData = obj.shapeLibrary.data(magShape).array;
517 compressed.num_samples = shapeData(1);
518 compressed.data = shapeData(2:end);
519 mag = mr.decompressShape(compressed);
520 shapeData = obj.shapeLibrary.data(phaseShape).array;
521 compressed.num_samples = shapeData(1);
522 compressed.data = shapeData(2:end);
523 phase = mr.decompressShape(compressed);
524 rf.signal = amplitude*mag.*exp(1j*2*pi*phase);
525 timeShape = libData(4);
526 if timeShape>0
527 shapeData = obj.shapeLibrary.data(timeShape).array;
528 compressed.num_samples = shapeData(1);
529 compressed.data = shapeData(2:end);
530 rf.t = mr.decompressShape(compressed)*obj.rfRasterTime;
531 rf.shape_dur=ceil((rf.t(end)-eps)/obj.rfRasterTime)*obj.rfRasterTime;
532 else
533 % generate default time raster on the fly
534 rf.t = ((1:length(rf.signal))-0.5)'*obj.rfRasterTime;
535 rf.shape_dur=length(rf.signal)*obj.rfRasterTime;
536 end
538 rf.center = libData(5); % new in v150
539 rf.delay = libData(6); % changed in v150
540 rf.freqPPM = libData(7); % changed in v150
541 rf.phasePPM = libData(8); % new changed v150
542 rf.freqOffset = libData(9); % changed in v150
543 rf.phaseOffset = libData(10); % new changed v150
545 rf.deadTime = obj.sys.rfDeadTime;
546 rf.ringdownTime = obj.sys.rfRingdownTime;
548% % SK: Is this a hack? (MZ: see below)
549% if length(libData) < 8
550% libData(8) = 0;
551% end
552% rf.deadTime = libData(9);
553% % SK: Using the same hack here
554% if length(libData) < 9
555% libData(9) = 0;
556% end
557% rf.ringdownTime = libData(9);
559 if nargin<=2
560 error('Parameter ''use'' is not optional since v1.5.0');
561 end
562 %TODO: fixme : use map built from mr.getSupportedRfUse();
563 switch use
564 case 'e'
565 rf.use='excitation';
566 case 'r'
567 rf.use='refocusing';
568 case 'i'
569 rf.use='inversion';
570 case 's'
571 rf.use='saturation';
572 case 'p'
573 rf.use='preparation';
574 otherwise
575 rf.use='undefined';
576 end
577 end
579 function [id shapeIDs]=registerRfEvent(obj, event)
580 % registerRfEvent : Add the event to the libraries (object,
581 % shapes, etc and return the event's ID. This ID should be
582 % stored in the object to accelerate addBlock()
584 mag = abs(event.signal);
585 amplitude = max(mag);
586 mag = mag / amplitude;
587 phase = angle(event.signal);
588 phase(phase < 0) = phase(phase < 0) + 2*pi;
589 phase = phase / (2*pi);
590 may_exist=true;
592 if isfield(event,'shapeIDs')
593 shapeIDs=event.shapeIDs;
594 else
595 shapeIDs=[0 0 0];
597 magShape = mr.compressShape(mag(:));
598 data = [magShape.num_samples magShape.data];
599 [shapeIDs(1),found] = obj.shapeLibrary.find_or_insert(data);
600 may_exist=may_exist & found;
602 phaseShape = mr.compressShape(phase);
603 data = [phaseShape.num_samples phaseShape.data];
604 [shapeIDs(2),found] = obj.shapeLibrary.find_or_insert(data);
605 may_exist=may_exist & found;
607 timeShape = mr.compressShape(event.t/obj.rfRasterTime); % time shape is stored in units of RF raster
608 if length(timeShape.data)==4 && all(timeShape.data == [0.5 1 1 timeShape.num_samples-3])
609 shapeIDs(3)=0;
610 else
611 data = [timeShape.num_samples timeShape.data];
612 [shapeIDs(3),found] = obj.shapeLibrary.find_or_insert(data);
613 may_exist=may_exist & found;
614 end
615 end
617 if isfield(event,'use')
618 % todo: fixme: use map from getSupportedRfUse
619 switch event.use
620 case {'excitation','refocusing','inversion','saturation','preparation','other'}
621 use = event.use(1);
622 otherwise
623 if strcmp(event.use,'u')
624 event.use='undefined'; % make it little more user-friendly
625 end
626 warning('Unknown or undefined RF pulse intended use ''use''=%s. Keep in mind that the ''use'' parameter is not optional since v1.5.0',event.use);
627 use = 'u'; % undefined
628 end
629 else
630 error('Parameter ''use'' is not optional since v1.5.0');
631 end
633 data = [amplitude shapeIDs(1) shapeIDs(2) shapeIDs(3) ...
634 event.center event.delay event.freqPPM event.phasePPM event.freqOffset event.phaseOffset ];%...
635 %event.deadTime event.ringdownTime];
636 if may_exist
637 id = obj.rfLibrary.find_or_insert(data,use);
638 else
639 id = obj.rfLibrary.insert(0,data,use);
640 end
642 if isfield(event,'name')
643 obj.rfID2NameMap(id) = event.name;
644 end
645 end
647 function [id,shapeIDs]=registerGradEvent(obj, event)
648 % registerGradEvent : Add the event to the libraries (object,
649 % shapes, etc and return the event's ID. This ID should be
650 % stored in the object to accelerate addBlock()
651 may_exist=true;
652 switch event.type
653 case 'grad'
654 amplitude = max(abs(event.waveform));
655 if amplitude>0
656 [~,~,fnz]=find(event.waveform,1); % find the first non-zero value and make it positive
657 amplitude=amplitude*sign(fnz);
658 end
659 if isfield(event,'shapeIDs')
660 shapeIDs=event.shapeIDs;
661 else
662 shapeIDs=[0 0];
663 % fill the shape IDs
664 if amplitude~=0
665 g = event.waveform./amplitude;
666 else
667 g = event.waveform;
668 end
669 c_shape = mr.compressShape(g);
670 s_data = [c_shape.num_samples c_shape.data];
671 [shapeIDs(1),found] = obj.shapeLibrary.find_or_insert(s_data);
672 may_exist=may_exist & found;
673 c_time = mr.compressShape(event.tt/obj.gradRasterTime);
674 if (length(c_time.data)==4 && all(c_time.data == [0.5 1 1 c_time.num_samples-3]))
675 % conventional grad on standard raster: shapeID
676 % is readily 0, nothing needs to be done
677 %shapeIDs(2)=0;
678 elseif (length(c_time.data)==3 && all(c_time.data == [0.5 0.5 c_time.num_samples-2]))
679 % grad on a half-raster (oversampling): shapeID
680 % needs to be set to -1 as a flag for oversampling
681 shapeIDs(2)=-1;
682 else
683 t_data = [c_time.num_samples c_time.data];
684 [shapeIDs(2),found] = obj.shapeLibrary.find_or_insert(t_data);
685 may_exist=may_exist & found;
686 end
687 end
688 data = [amplitude event.first event.last shapeIDs event.delay];
689 case 'trap'
690 data = [event.amplitude event.riseTime ...
691 event.flatTime event.fallTime ...
692 event.delay];
693 otherwise
694 error('unknown grdient type passed to registerGradEvent()');
695 end
696 if may_exist
697 id = obj.gradLibrary.find_or_insert(data,event.type(1));
698 else
699 id = obj.gradLibrary.insert(0,data,event.type(1));
700 end
702 if isfield(event,'name')
703 obj.gradID2NameMap(id) = event.name;
704 end
705 end
707 function [id,shapeID]=registerAdcEvent(obj, event)
708 % registerAdcEvent : Add the event to the libraries (object,
709 % shapes, etc and return the event's ID. This ID should be
710 % stored in the object to accelerate addBlock()
712 surely_new=false;
713 if isempty(event.phaseModulation)
714 shapeID=0;
715 else
716 if isfield(event,'shapeID')
717 shapeID=event.shapeID;
718 else
719 phaseShape = mr.compressShape(event.phaseModulation(:));
720 data = [phaseShape.num_samples phaseShape.data];
721 [shapeID,shape_found] = obj.shapeLibrary.find_or_insert(data);
722 if ~shape_found
723 surely_new=true;
724 end
725 end
726 end
728 data = [event.numSamples event.dwell max(event.delay,event.deadTime), ... % MZ: replaced event.delay+event.deadTime with a max(...) because we allow for overlap of the delay and the dead time
729 event.freqPPM event.phasePPM event.freqOffset event.phaseOffset shapeID]; % event.deadTime];
730 if surely_new
731 id = obj.adcLibrary.insert(0,data);
732 else
733 id = obj.adcLibrary.find_or_insert(data);
734 end
736 if isfield(event,'name')
737 obj.adcID2NameMap(id) = event.name;
738 end
739 end
741 function id=registerControlEvent(obj, event)
742 % registerControlEvent : Add the event to the libraries (object,
743 % shapes, etc and return the event's ID. This ID should be
744 % stored in the object to accelerate addBlock()
745 event_type=find(strcmp(event.type,{'output','trigger'}));
746 if (event_type==1)
747 event_channel=find(strcmp(event.channel,{'osc0','osc1','ext1'})); % trigger codes supported by the Siemens interpreter as of May 2019
748 elseif (event_type==2)
749 event_channel=find(strcmp(event.channel,{'physio1','physio2'})); % trigger codes supported by the Siemens interpreter as of June 2019
750 else
751 error('unsupported control event type');
752 end
753 data = [event_type event_channel event.delay event.duration];
754 id = obj.trigLibrary.find_or_insert(data);
755 end
757 function id=registerLabelEvent(obj, event)
758 % registerLabelEvent : Add the event to the libraries (object,
759 % shapes, etc and return the event's ID. This ID should be
760 % stored in the object to accelerate addBlock()
761 label_id=find(strcmp(event.label,mr.getSupportedLabels()));
762 data=[event.value label_id];
763 switch event.type
764 case 'labelset'
765 id = obj.labelsetLibrary.find_or_insert(data);
766 case 'labelinc'
767 id = obj.labelincLibrary.find_or_insert(data);
768 otherwise
769 error('unknown label type passed to registerLabelEvent()');
770 end
771 end
773 function id=registerSoftDelayEvent(obj, event)
774 % registerDeleyEvent : Add the event to the libraries (object,
775 % shapes, etc and return the event's ID. This ID should be
776 % stored in the object to accelerate addBlock()
777 try
778 hintID=obj.softDelayHints1(event.hint);
779 catch
780 hintID=obj.softDelayHints1.length()+1;
781 obj.softDelayHints1(event.hint)=hintID;
782 obj.softDelayHints2{hintID}=event.hint;
783 end
784 data = [event.num event.offset event.factor hintID];
785 id = obj.softDelayLibrary.find_or_insert(data);
786 end
788 function id=registerRfShimEvent(obj, event)
789 % registerRfShimEvent : Add the event to the libraries (object,
790 % shapes, etc and return the event's ID. This ID should be
791 % stored in the object to accelerate addBlock()
792 data = [abs(event.shimVector(:)),angle(event.shimVector(:))].'; % make data(:) to produce an interleaving vector of amplitudes and phases
793 id = obj.rfShimLibrary.find_or_insert(data(:));
794 end
796 function id=registerRotationEvent(obj, event)
797 % registerRotationEvent : Add the event to the libraries (object,
798 % shapes, etc and return the event's ID. This ID should be
799 % stored in the object to accelerate addBlock()
800 data = event.rotQuaternion;
801 % confirm that the rotation matrix is valid
802 if ( length(data) == 4 ) & ( abs(1.0 - sum(data.^2)) < 1e-6 )
803 id = obj.rotationLibrary.find_or_insert(data(:)');
804 else
805 error('invalid rotation quaternion detected during registerRotationEvent()');
806 end
807 end
809 %TODO: Replacing blocks in the middle of sequence can cause unused
810 %events in the libraries. These can be detected and pruned.
811 function setBlock(obj, index, varargin)
812 %setBlock Replace or add sequence block.
813 % setBlock(obj, index, bStruct) Replace block at index with new
814 % block provided as block structure.
815 %
816 % setBlock(obj, index, e1, e2, ...) Create a new block from
817 % events and store at position given by index.
818 %
819 % setBlock(obj, index, duration, e1, e2, ...) Create a new
820 % block with the given predefined duration populated with
821 % events e1, e2, etc. and store at position given by index.
822 % If the duration of any of the events exceeds the desired
823 % duration an error will be thrown.
824 %
825 % The block or events are provided in uncompressed form and
826 % will be stored in the compressed, non-redundant internal
827 % libraries.
828 %
829 % See also getBlock, addBlock
831 % Convert block structure to cell array of events
832 varargin=mr.block2events(varargin);
834 newBlock=zeros(1,7);
835 duration = 0;
837 check_g = cell(1,3); % cell-array containing a structure, each with the index and pairs of gradients/times
838 extensions = [];
839 required_duration=[];
840 rotQuaternion = []; % rotation extension
841 roundUpBlockDuration=false; % historical default
843 % Loop over events adding to library if necessary and creating
844 % block event structure.
845 for i = 1:length(varargin)
846 event = varargin{i};
847 if isstruct(event)
848 switch event(1).type % we accept multiple extensions and one of the possibilities is an array of extensions
849 case 'rf'
850 if isfield(event,'id')
851 newBlock(2)=event.id;
852 else
853 newBlock(2) = obj.registerRfEvent(event);
854 end
855 duration = max(duration, event.shape_dur + event.delay + event.ringdownTime);
856 case 'grad'
857 channelNum = find(strcmp(event.channel, ...
858 {'x', 'y', 'z'}));
860 idx = 2 + channelNum;
861 grad_duration = event.delay + ceil(event.tt(end)/obj.gradRasterTime-1e-10)*obj.gradRasterTime;
863 grad_start = event.delay + floor(event.tt(1)/obj.gradRasterTime+1e-10)*obj.gradRasterTime;
865 %check_g{channelNum}.idx = idx;
866 check_g{channelNum}.start = [grad_start, event.first];
867 check_g{channelNum}.stop = [grad_duration, event.last];
869 if newBlock(idx)>0
870 error('Trying to add more than one gradient per axis on axis %s in block %d',event.channel,index);
871 end
872 if isfield(event,'id')
873 newBlock(idx) = event.id;
874 else
875 newBlock(idx) = obj.registerGradEvent(event);
876 end
877 duration = max(duration, grad_duration);
879 case 'trap'
880 channelNum = find(strcmp(event.channel,{'x','y','z'}));
883 idx = 2 + channelNum;
885 % MZ: the checks as implemented below only make sense for non-trapezoid gradients, commented out the coe below
886 % %check_g{channelNum}.idx = idx;
887 % check_g{channelNum}.start = [0, 0];
888 % check_g{channelNum}.stop = [event.delay + ...
889 % event.riseTime + ...
890 % event.fallTime + ...
891 % event.flatTime, 0];
893 if newBlock(idx)>0
894 error('Trying to add more than one gradient per axis on axis %s in block %d',event.channel,index);
895 end
896 if isfield(event,'id')
897 newBlock(idx) = event.id;
898 else
899 newBlock(idx) = obj.registerGradEvent(event);
900 end
901 duration=max(duration,event.delay+event.riseTime+event.flatTime+event.fallTime);
903 case 'adc'
904 if isfield(event,'id')
905 newBlock(6) = event.id;
906 else
907 newBlock(6) = obj.registerAdcEvent(event);
908 end
909 duration=max(duration,event.delay+event.numSamples*event.dwell+event.deadTime); % adcDeadTime is added after the sampling period (mr.makeADC also adds a delay before the actual sampling if it was shorter)
910 case 'delay'
911 %if isfield(event,'id')
912 % id=event.id;
913 %else
914 % id = obj.registerDelayEvent(event);
915 %end
916 %newBlock(1)=id;
917 % delay is not a true event any more so we account
918 % for the duration but do not add anything
919 duration=max(duration,event.delay);
920 case {'output','trigger'}
921 for e=event % allow multiple extensions as an array
922 if isfield(e,'id')
923 id=e.id;
924 else
925 id=obj.registerControlEvent(e);
926 end
927 %newBlock(7)=id; % now we just
928 % collect the list of extension objects and we will
929 % add it to the event table later
930 % ext=struct('type', 1, 'ref', id);
931 ext=struct('type', obj.getExtensionTypeID('TRIGGERS'), 'ref', id);
932 extensions=[extensions ext];
933 duration=max(duration,e.delay+e.duration);
934 end
935 case {'labelset','labelinc'}
936 for e=event % allow multiple extensions as an array
937 if isfield(e,'id')
938 id=e.id;
939 else
940 id=obj.registerLabelEvent(e);
941 end
942 % % label_id=find(strcmp(e.label,mr.getSupportedLabels()));
943 % % data=[e.value label_id];
944 % % [id,found] = obj.labelsetLibrary.find(data);
945 % % if ~found
946 % % obj.labelsetLibrary.insert(id,data);
947 % % end
949 % collect the list of extension objects and we will
950 % add it to the event table later
951 %ext=struct('type', 2, 'ref', id);
952 ext=struct('type', obj.getExtensionTypeID(upper(e.type)), 'ref', id);
953 extensions=[extensions ext];
954 end
955 case 'softDelay'
956 for e=event % allow multiple extensions as an array
957 if isfield(e,'id')
958 id=e.id;
959 else
960 id=obj.registerSoftDelayEvent(e);
961 end
962 % collect the list of extension objects and
963 % add it to the event table later
964 % ext=struct('type', 1, 'ref', id);
965 ext=struct('type', obj.getExtensionTypeID('DELAYS'), 'ref', id);
966 extensions=[extensions ext];
967 end
968 case 'rfShim'
969 if isfield(event,'id') % this is an array if we have a loop?
970 id=event.id;
971 else
972 id=obj.registerRfShimEvent(event);
973 end
974 % collect the list of extension objects and
975 % add it to the event table later
976 % ext=struct('type', 1, 'ref', id);
977 ext=struct('type', obj.getExtensionTypeID('RF_SHIMS'), 'ref', id);
978 extensions=[extensions ext];
979 case 'rot3D'
980 if ~isempty(rotQuaternion)
981 error('Only one ''rotation'' extension event can be added per block');
982 end
983 rotQuaternion=event.rotQuaternion;
984 if isfield(event,'id')
985 id=event.id;
986 else
987 id = obj.registerRotationEvent(event);
988 end
989 ext=struct('type', obj.getExtensionTypeID('ROTATIONS'), 'ref', id);
990 extensions=[extensions ext];
991 otherwise
992 error('Attempting to add an unknown event to the block.');
993 end
994 else
995 if isnumeric(event)
996 % interpret the single numeric parameter as a
997 % requested duration, but throw an error if
998 % multiple numbers are provided
999 if isempty(required_duration)
1000 required_duration=event;
1001 else
1002 error('More than one numeric parameter given to setBlock()');
1003 end
1004 elseif ischar(event) && strcmp(event,'roundUpBlockDuration')
1005 roundUpBlockDuration=true;
1006 else
1007 warning('Unknown parameter passed to block %d', index);
1008 end
1009 end
1010 end
1012 if ~isempty(extensions)
1013 % add extensions now... but it's tricky actually
1014 % we need to check whether the exactly the same list if
1015 % extensions already exists, otherwise we have to create a
1016 % new one... ooops, we have a potential problem with the
1017 % key mapping then... The trick is that we rely on the
1018 % sorting of the extension IDs and then we can always find
1019 % the last one in the list by setting the reference to the
1020 % next to 0 and then proceed with the other elements.
1021 [~,I]=sort([extensions(:).ref]);
1022 extensions=extensions(I);
1023 all_found=true;
1024 id=0;
1025 for i=1:length(extensions)
1026 data=[extensions(i).type extensions(i).ref id];
1027 [id,found] = obj.extensionLibrary.find(data);
1028 all_found = all_found && found;
1029 if ~found
1030 break;
1031 end
1032 end
1033 if ~all_found
1034 % add the list
1035 id=0;
1036 for i=1:length(extensions)
1037 data=[extensions(i).type extensions(i).ref id];
1038 [id,found] = obj.extensionLibrary.find(data);
1039 if ~found
1040 obj.extensionLibrary.insert(id,data);
1041 end
1042 end
1043 end
1044 % sanity checks for the softDelay
1045 nSoftDelays=sum([extensions(:).type]==obj.getExtensionTypeID('DELAYS'));
1046 if nSoftDelays
1047 if nSoftDelays>1
1048 error('Only one ''softDelay'' extension event can be added per block');
1049 end
1050 if duration==0 && isempty(required_duration)
1051 error('Soft delay extension can only be used in conjunstion with blocks of non-zero duration'); % otherwise the gradient checks get tedious
1052 end
1053 if any(newBlock(2:6)~=0)
1054 error('Soft delay extension can only be used in empty blocks (blocks containing no conventional events such as RF, adc or gradients).')
1055 end
1056 end
1057 % now we add the ID
1058 newBlock(7)=id;
1059 end
1061 if duration>0
1063 if roundUpBlockDuration
1064 duration=ceil(duration/obj.sys.blockDurationRaster)*obj.sys.blockDurationRaster;
1065 end
1067 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1068 %%% PERFORM GRADIENT CHECKS %%%
1069 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1071 % see if we have the valid preceding gradient value data
1072 %gradCheckData % struct('validForBlockNum',0,'lastGradVals', [0 0 0]);
1074 if index > 1 && obj.gradCheckData.validForBlockNum ~= index-1
1075 % need to update gradCheckData
1076 obj.gradCheckData.validForBlockNum = index-1;
1077 obj.gradCheckData.lastGradVals(:)=0;
1078 [~,prev_nonempty_block]=find(obj.blockDurations(1:(index-1))>0, 1, 'last');
1079 if ~isempty(prev_nonempty_block)
1080 for i= 1:length(obj.gradCheckData.lastGradVals) % TODO: MZ: check this with external gradient channels !!!
1081 prev_id = obj.blockEvents{prev_nonempty_block}(i+2); % careful! direct eventLib access
1082 if prev_id ~= 0
1083 prev_lib = obj.gradLibrary.get(prev_id); % MZ: for performance reasons we access the gradient library directly. I know, this is not elegant
1084 prev_dat = prev_lib.data;
1085 prev_type = prev_lib.type;
1086 if prev_type == 'g'
1087 obj.gradCheckData.lastGradVals(i) = prev_dat(3); % change in v150 - now it is '3' % '6' means last; MZ: I know, this is a real hack...
1088 end
1089 end
1090 end
1091 end
1092 % now it's a pain, but we also need to extract the rotation extension of the previous block --
1093 warning('mr:fixmePreviousRotationExtension','FIXME: need rotation extension of the previous non-zero block here...');
1094 end
1095 % check if connection to the previous block is correct using check_g and gradCheckData.lastGradVals
1096 % up to here gradCheckData.lastGradVals are in physical coordinates, we transform them into current
1097 % logical coordinates of this block has rotation extension
1098 if ~isempty(rotQuaternion)
1099 obj.gradCheckData.lastGradVals = mr.aux.quat.rotate(mr.aux.quat.conjugate(rotQuaternion), obj.gradCheckData.lastGradVals);
1100 end
1101 for i= 1:3 %length(check_g) % TODO: MZ: check this with external gradient channels !!!
1102 cg=check_g{i}; % cg_temp is still a cell-array with a single element here...
1103 % connection to the previous block in case of extended or shaped gradients
1104 if isempty(cg)
1105 if abs(obj.gradCheckData.lastGradVals(i)) > obj.sys.maxSlew * obj.sys.gradRasterTime
1106 error('Error in block %d on gradient axis %d: previous block ended with non-zero amplitude but the current block has no compatible gradient.', index, i);
1107 end
1108 % update the gradCheckData.lastGradVals(i)
1109 obj.gradCheckData.lastGradVals(i)=0;
1110 continue;
1111 end
1113 % check the start
1114 if abs(cg.start(2)) > obj.sys.maxSlew * obj.sys.gradRasterTime % MZ: we only need the following check if the current gradient starts at non-0
1115 if cg.start(1) ~= 0
1116 error('Error in block %d: No delay allowed for gradients which start with a non-zero amplitude', index);
1117 end
1118 if index > 1
1119 if abs(obj.gradCheckData.lastGradVals(i) - cg.start(2)) > obj.sys.maxSlew * obj.sys.gradRasterTime
1120 error('Error in block %d on gradient axis %d: Two consecutive gradients need to have the same amplitude at the connection point', index, i);
1121 end
1122 else
1123 error('First gradient in the the first block has to start at 0.');
1124 end
1125 end
1127 % Check if gradients, which do not end at 0, are as long as the block itself.
1128 if cg.stop(2) > obj.sys.maxSlew * obj.sys.gradRasterTime && abs(cg.stop(1)-duration) > 1e-7
1129 error('Error in block %d: A gradient that doesn''t end at zero needs to be aligned to the block boundary', index);
1130 end
1132 % update the gradCheckData.lastGradVals(i)
1133 obj.gradCheckData.lastGradVals(i)=cg.stop(2); % we can play with the continuity check values by transforming them back and foth, here to physical, at the baginning of the check to current logical
1134 end
1135 % now transfrom the gradCheckData.lastGradVals to physical coordinates if the present block has rotation
1136 if ~isempty(rotQuaternion)
1137 obj.gradCheckData.lastGradVals = mr.aux.quat.rotate(rotQuaternion, obj.gradCheckData.lastGradVals);
1138 end
1139 end
1140 % finish updating gradCheckData (if current block duration is 0 we simply update the validity indicator)
1141 obj.gradCheckData.validForBlockNum = index;
1143 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1144 %%% GRADIENT CHECKS DONE %%%
1145 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1148 % now copy the block into the internal data structure
1149 obj.blockEvents{index}=newBlock;
1151 if ~isempty(required_duration)
1152 if duration-required_duration>eps
1153 error('Required block duration is %g s but actuall block duration is %g s', required_duration, duration);
1154 end
1155 duration=required_duration;
1156 end
1158 obj.blockDurations(index)=duration;
1159 end
1161 function raw_block = getRawBlockContentIDs(obj, index)
1162 %getRawBlockContentIDs Return a block content of the sequence.
1163 % b=getRawBlockContentIDs(obj, index) Return the block
1164 % content IDs specified by the index of the block.
1165 %
1166 % No block events are created, only the IDs of the objects
1167 % are returned.
1168 %
1169 % See also getBlock, setBlock, addBlock
1171 raw_block=struct('blockDuration', 0, 'rf', [], 'gx', [], 'gy', [], 'gz', [], 'adc', [], 'ext', [] );
1173 eventInd = obj.blockEvents{index};
1175 if eventInd(7) > 0
1176 % we have extensions -- triggers, labels, etc
1177 % first find how many and preallocate raw_block.ext array
1178 nextExtID=eventInd(7);
1179 cExt=0;
1180 while nextExtID~=0
1181 cExt=cExt+1;
1182 % now update nextExtID
1183 nextExtID=obj.extensionLibrary.data(nextExtID).array(3);
1184 end
1185 raw_block.ext=zeros(2,cExt);
1186 % now scan through the list again and extract the extension data
1187 nextExtID=eventInd(7);
1188 cExt=0;
1189 while nextExtID~=0
1190 cExt=cExt+1;
1191 extData = obj.extensionLibrary.data(nextExtID).array;
1192 % format: extType, extID, nextExtID
1193 raw_block.ext(:,cExt)=extData(1:2);
1194 % now update nextExtID
1195 nextExtID=extData(3);
1196 end
1197 end
1198 if eventInd(2) > 0
1199 raw_block.rf=eventInd(2);
1200 end
1201 gradChannels = {'gx', 'gy', 'gz'};
1202 for i = 1:length(gradChannels)
1203 if eventInd(2+i) > 0
1204 raw_block.(gradChannels{i})=eventInd(2+i);
1205 end
1206 end
1207 if eventInd(6) > 0
1208 raw_block.adc = eventInd(6);
1209 end
1210 end
1212 function block = getBlock(obj, index, addIDs)
1213 %getBlock Return a block of the sequence.
1214 % b=getBlock(obj, index) Return the block specified by the
1215 % index.
1216 %
1217 % The block is created from the sequence data with all
1218 % events and shapes decompressed.
1219 %
1220 % See also setBlock, addBlock
1222 if nargin < 3
1223 addIDs=false;
1224 end
1226 block=struct('blockDuration', 0, 'rf', [], 'gx', [], 'gy', [], 'gz', [], 'adc', [] );
1228 %block(1).rf = [];
1229 %eventInd = obj.blockEvents{index};
1230 raw_block = obj.getRawBlockContentIDs(index);
1232 if ~isempty(raw_block.ext)
1233 % we have extensions -- triggers, labels, etc
1234 % ext field format: extType, extID
1236 % unpack trigger(s)
1237 trig_ext=raw_block.ext(2,raw_block.ext(1,:)==obj.getExtensionTypeID('TRIGGERS'));
1238 if ~isempty(trig_ext)
1239 trigger_types={'output','trigger'};
1240 for i=length(trig_ext):-1:1 % backwards for preallocation
1241 data = obj.trigLibrary.data(trig_ext(i)).array;
1242 trig.type = trigger_types{data(1)};
1243 if (data(1)==1)
1244 trigger_channels={'osc0','osc1','ext1'};
1245 trig.channel=trigger_channels{data(2)};
1246 elseif (data(1)==2)
1247 trigger_channels={'physio1','physio2'};
1248 trig.channel=trigger_channels{data(2)};;
1249 else
1250 error('unsupported trigger event type');
1251 end
1252 trig.delay = data(3);
1253 trig.duration = data(4);
1254 if addIDs
1255 trig.id=trig_ext(i);
1256 end
1257 % we allow for multiple triggers per block
1258 block.trig(i)=trig;
1259 end
1260 end
1261 % unpack labels
1262 lid_set=obj.getExtensionTypeID('LABELSET');
1263 lid_inc=obj.getExtensionTypeID('LABELINC');
1264 supported_labels=mr.getSupportedLabels();
1265 label_ext=raw_block.ext(:,raw_block.ext(1,:)==lid_set | raw_block.ext(1,:)==lid_inc);
1266 if ~isempty(label_ext)
1267 for i=size(label_ext,2):-1:1 % backwards for preallocation
1268 if label_ext(1,i)==lid_set
1269 label.type='labelset';
1270 data = obj.labelsetLibrary.data(label_ext(2,i)).array;
1271 else
1272 label.type='labelinc';
1273 data = obj.labelincLibrary.data(label_ext(2,i)).array;
1274 end
1275 label.label=supported_labels{data(2)};
1276 label.value=data(1);
1277 if addIDs
1278 label.id=label_ext(2,i);
1279 end
1280 block.label(i) = label;
1281 end
1282 end
1283 % unpack RF shim
1284 rf_shim_ext=raw_block.ext(:,raw_block.ext(1,:)==obj.getExtensionTypeID('RF_SHIMS'));
1285 if ~isempty(rf_shim_ext)
1286 if size(rf_shim_ext,2)>1
1287 error('Only one RF shim extension object per block is allowed');
1288 end
1289 data = obj.rfShimLibrary.data(rf_shim_ext(2,1)).array;
1290 if addIDs
1291 block.rfShim=struct('type','rfShim','shimVector',data(1:2:end).*exp(1i*data(2:2:end)),'id',rf_shim_ext(2,1));
1292 else
1293 block.rfShim=struct('type','rfShim','shimVector',data(1:2:end).*exp(1i*data(2:2:end)));
1294 end
1295 end
1296 % unpack 3D rotations
1297 rotation_ext=raw_block.ext(:,raw_block.ext(1,:)==obj.getExtensionTypeID('ROTATIONS'));
1298 if ~isempty(rotation_ext)
1299 if size(rotation_ext,2)>1
1300 error('Only one rotation extension object per block is allowed');
1301 end
1302 data = obj.rotationLibrary.data(rotation_ext(2,1)).array;
1303 if addIDs
1304 block.rotation=struct('type','rot3D','rotQuaternion',data,'id',rotation_ext(2,1));
1305 else
1306 block.rotation=struct('type','rot3D','rotQuaternion',data);
1307 end
1308 end
1309 % unpack delay
1310 delay_ext=raw_block.ext(:,raw_block.ext(1,:)==obj.getExtensionTypeID('DELAYS'));
1311 if ~isempty(delay_ext)
1312 if size(delay_ext,2)>1
1313 error('Only one soft delay extension object per block is allowed');
1314 end
1315 data = obj.softDelayLibrary.data(delay_ext(2,1)).array;
1316 if addIDs
1317 block.softDelay=struct('type','softDelay','num',data(1),'offset',data(2),'factor',data(3),'hint',obj.softDelayHints2{data(4)},'id',delay_ext(2,1));
1318 else
1319 block.softDelay=struct('type','softDelay','num',data(1),'offset',data(2),'factor',data(3),'hint',obj.softDelayHints2{data(4)});
1320 end
1321 end
1322 if length(trig_ext)+size(label_ext,2)~=size(raw_block.ext,2)
1323 for i=1:size(raw_block.ext,2)
1324 if raw_block.ext(1,i)~=obj.getExtensionTypeID('TRIGGERS') && ...
1325 raw_block.ext(1,i)~=obj.getExtensionTypeID('LABELSET') && ...
1326 raw_block.ext(1,i)~=obj.getExtensionTypeID('LABELSET') && ...
1327 raw_block.ext(1,i)~=obj.getExtensionTypeID('RF_SHIMS') && ...
1328 raw_block.ext(1,i)~=obj.getExtensionTypeID('DELAYS') && ...
1329 raw_block.ext(1,i)~=obj.getExtensionTypeID('ROTATIONS')
1330 warning('unknown extension ID %d', raw_block.ext(1,i));
1331 end
1332 end
1333 end
1334 end
1335 % finished extensions
1337 if ~isempty(raw_block.rf)
1338 if length(obj.rfLibrary.type)>=raw_block.rf
1339 block.rf = obj.rfFromLibData(obj.rfLibrary.data(raw_block.rf).array,obj.rfLibrary.type(raw_block.rf));
1340 else
1341 block.rf = obj.rfFromLibData(obj.rfLibrary.data(raw_block.rf).array); % undefined type/use
1342 end
1343 if addIDs
1344 block.rf.id=raw_block.rf;
1345 % this is a bit of a hack because we have to access the rfLibrary directly
1346 block.rf.shapeIDs=obj.rfLibrary.data(raw_block.rf).array(2:4); % ampl_shape_id phase_shape_id time_shape_id
1347 end
1348 end
1349 gradChannels = {'gx', 'gy', 'gz'};
1350 for i = 1:length(gradChannels)
1351 gid=raw_block.(gradChannels{i});
1352 if ~isempty(gid)
1353 type = obj.gradLibrary.type(gid);
1354 libData = obj.gradLibrary.data(gid).array;
1355 grad=struct();
1356 if type == 't'
1357 grad.type = 'trap';
1358 else
1359 grad.type = 'grad';
1360 end
1361 grad.channel = gradChannels{i}(2);
1362 if strcmp(grad.type,'grad')
1363 amplitude = libData(1);
1364 shapeId = libData(4); % change in v150
1365 timeId = libData(5); % change in v150
1366 delay = libData(6); % change in v150
1367 shapeData = obj.shapeLibrary.data(shapeId).array;
1368 compressed.num_samples = shapeData(1);
1369 compressed.data = shapeData(2:end);
1370 try
1371 g = mr.decompressShape(compressed);
1372 catch
1373 fprintf(' mr.decompressShape() failed for shapeId %d\n', shapeId);
1374 error('mr.decompressShape() failed for shapeId %d', shapeId);
1375 end
1376 grad.waveform = amplitude*g;
1377 % SK: This looks like a bug to me.
1378% grad.t = (1:length(g))'*obj.gradRasterTime;
1379 if (timeId==0)
1380 grad.tt = ((1:length(g))-0.5)'*obj.gradRasterTime; % TODO: evetually we may remove these true-times
1381 t_end=length(g)*obj.gradRasterTime;
1382 %grad.t = (0:length(g)-1)'*obj.gradRasterTime;
1383 grad.area=sum(grad.waveform)*obj.gradRasterTime;
1384 elseif (timeId==-1)
1385 % gradient with oversampling by a factor of 2
1386 grad.tt = ((1:length(g)))'/2*obj.gradRasterTime;
1387 assert(length(grad.tt)==length(grad.waveform));
1388 assert(mod(length(g),2)==1);
1389 t_end=(length(g)+1)/2*obj.gradRasterTime;
1390 grad.area=sum(grad.waveform(1:2:end))*obj.gradRasterTime; % remove oversampling
1391 else
1392 tShapeData = obj.shapeLibrary.data(timeId).array;
1393 compressed.num_samples = tShapeData(1);
1394 compressed.data = tShapeData(2:end);
1395 try
1396 grad.tt = mr.decompressShape(compressed)*obj.gradRasterTime;
1397 catch
1398 fprintf(' mr.decompressShape() failed for shapeId %d\n', shapeId);
1399 error('mr.decompressShape() failed for shapeId %d', shapeId);
1400 end
1401 assert(length(grad.waveform) == length(grad.tt));
1402 t_end=grad.tt(end);
1403 grad.area=0.5*sum((grad.tt(2:end)-grad.tt(1:end-1)).*(grad.waveform(2:end)+grad.waveform(1:end-1)));
1404 end
1405 grad.shape_id=shapeId; % needed for the second pass of read()
1406 grad.time_id=timeId; % needed for the second pass of read()
1407 grad.delay = delay;
1408 grad.shape_dur = t_end;
1409 grad.first = libData(2); % change in v150 - we always have first/last now
1410 grad.last = libData(3); % change in v150 - we always have first/last now
1411 if addIDs
1412 grad.shapeIDs = [shapeId timeId];
1413 end
1414 else
1415 grad.amplitude = libData(1);
1416 grad.riseTime = libData(2);
1417 grad.flatTime = libData(3);
1418 grad.fallTime = libData(4);
1419 grad.delay = libData(5);
1420 grad.area = grad.amplitude*(grad.flatTime + ...
1421 grad.riseTime/2 + ...
1422 grad.fallTime/2);
1423 grad.flatArea = grad.amplitude*grad.flatTime;
1424 end
1425 if addIDs
1426 grad.id=gid;
1427 end
1428 block.(gradChannels{i}) = grad;
1429 end
1430 end
1431 if ~isempty(raw_block.adc)
1432 libData = obj.adcLibrary.data(raw_block.adc).array;
1433 shapeIdPhaseModulation=libData(end);
1434 if shapeIdPhaseModulation
1435 shapeData = obj.shapeLibrary.data(shapeIdPhaseModulation).array;
1436 compressed.num_samples = shapeData(1);
1437 compressed.data = shapeData(2:end);
1438 try
1439 phaseShape = mr.decompressShape(compressed);
1440 catch
1441 error('mr.decompressShape() failed for shapeId %d', shapeIdPhaseModulation);
1442 end
1443 else
1444 phaseShape=0; % wee need a 0 trick here because [] did not work
1445 end
1446 adc = cell2struct(num2cell([libData(1:end-1) 0 obj.sys.adcDeadTime]), ...
1447 {'numSamples', 'dwell', 'delay', ...
1448 'freqPPM', 'phasePPM', 'freqOffset', 'phaseOffset', ...
1449 'phaseModulation','deadTime'}, 2);
1450 if shapeIdPhaseModulation
1451 adc.phaseModulation=phaseShape;
1452 else
1453 adc.phaseModulation=[]; % replace 0 with an empty array
1454 end
1455 adc.type = 'adc';
1456 if addIDs
1457 adc.id=raw_block.adc;
1458 end
1459 block.adc = adc;
1460 end
1461 block.blockDuration=obj.blockDurations(index);
1462% % now that delays in v1.4 and later Pulseq revisions are not
1463% % stored, we need to see whether we need a delay to explain the
1464% % current block duration
1465% if (mr.calcDuration(block)~=obj.blockDurations(index))
1466% tmpDelay.type = 'delay';
1467% tmpDelay.delay = obj.blockDurations(index);
1468% block.delay = tmpDelay;
1469% end
1470 end
1472 function [ktraj_adc, ktraj, t_excitation, t_refocusing, t_adc] = calculateKspaceUnfunc(obj, varargin)
1473 % calculate the k-space trajectory of the entire pulse sequence
1474 % optional parameter 'trajectory_delay' sets the compensation
1475 % factor to align ADC and gradients in the reconstruction
1476 % Return values: ktraj_adc, ktraj, t_excitation, t_refocusing
1478 persistent parser
1479 if isempty(parser)
1480 parser = inputParser;
1481 parser.FunctionName = 'calculateKspace';
1482 parser.addParamValue('trajectory_delay',0,@(x)(isnumeric(x)));
1483 end
1484 parse(parser,varargin{:});
1485 opt = parser.Results;
1487 if any(abs(opt.trajectory_delay)>100e-6)
1488 warning('trajectory delay of (%s) us is suspiciously high',num2str(opt.trajectory_delay*1e6));
1489 end
1491 % initialise the counters and accumulator objects
1492 c_excitation=0;
1493 c_refocusing=0;
1494 c_adcSamples=0;
1495 % loop throught the blocks to prepare preallocations
1496 for iB=1:length(obj.blockEvents)
1497 block = obj.getBlock(iB);
1498 if ~isempty(block.rf)
1499 if (~isfield(block.rf,'use') || strcmp(block.rf.use,'excitation') || strcmp(block.rf.use,'undefined'))
1500 c_excitation=c_excitation+1;
1501 elseif strcmp(block.rf.use,'refocusing')
1502 c_refocusing=c_refocusing+1;
1503 end
1504 end
1505 if ~isempty(block.adc)
1506 c_adcSamples=c_adcSamples+block.adc.numSamples;
1507 end
1508 end
1510 %
1511 t_excitation=zeros(c_excitation,1);
1512 t_refocusing=zeros(c_refocusing,1);
1513 ktime=zeros(c_adcSamples,1);
1514 current_dur=0;
1515 c_excitation=1;
1516 c_refocusing=1;
1517 kcouter=1;
1518 traj_recon_delay=opt.trajectory_delay;
1520 % go through the blocks and collect RF and ADC timing data
1521 for iB=1:length(obj.blockEvents)
1522 block = obj.getBlock(iB);
1523 if ~isempty(block.rf)
1524 rf=block.rf;
1525 t=rf.delay+mr.calcRfCenter(rf);
1526 if (~isfield(block.rf,'use') || strcmp(block.rf.use,'excitation') || strcmp(block.rf.use,'undefined'))
1527 t_excitation(c_excitation) = current_dur+t;
1528 c_excitation=c_excitation+1;
1529 elseif strcmp(block.rf.use,'refocusing')
1530 t_refocusing(c_refocusing) = current_dur+t;
1531 c_refocusing=c_refocusing+1;
1532 end
1533 end
1534 if ~isempty(block.adc)
1535 ktime(kcouter:(kcouter-1+block.adc.numSamples)) = ((0:(block.adc.numSamples-1))+0.5)... % according to the information from Klaus Scheffler and indirectly from Siemens this is the present convention (the samples are shifted by 0.5 dwell)
1536 *block.adc.dwell + block.adc.delay + current_dur + traj_recon_delay;
1537 kcouter=kcouter+block.adc.numSamples;
1538 end
1539 current_dur=current_dur+obj.blockDurations(iB);%mr.calcDuration(block);
1540 end
1542 % now calculate the actual k-space trajectory based on the
1543 % gradient waveforms
1544 gw=obj.gradient_waveforms();
1545 i_excitation=round(t_excitation/obj.gradRasterTime);
1546 i_refocusing=round(t_refocusing/obj.gradRasterTime);
1547% ii_next_excitation=min(length(i_excitation),1);
1548% ii_next_refocusing=min(length(i_refocusing),1);
1549% ktraj=zeros(size(gw));
1550% k=[0;0;0];
1551% % TODO: replace this plain stupid loop with a segment-wise
1552% % integration (with segments defined by the RF pulses)
1553% for i=1:size(gw,2)
1554% k=k+gw(:,i)*obj.gradRasterTime;
1555% ktraj(:,i)=k;
1556% %if find(i_excitation==i,1)
1557% if ii_next_excitation>0 && i_excitation(ii_next_excitation)==i
1558% k=0;
1559% ktraj(:,i)=NaN; % we use NaN-s to mark the excitation point, they interrupt the plots
1560% ii_next_excitation = min(length(i_excitation),ii_next_excitation+1);
1561% end
1562% %if find(i_refocusing==i,1)
1563% if ii_next_refocusing>0 && i_refocusing(ii_next_refocusing)==i
1564% k=-k;
1565% ii_next_refocusing = min(length(i_refocusing),ii_next_refocusing+1);
1566% end
1567% end
1568 i_periods=sort([1; i_excitation+1; i_refocusing+1; size(gw,2)+1]); % we need thise +1 for compatibility with the above code which prooved to be correct
1569 ii_next_excitation=min(length(i_excitation),1);
1570 ii_next_refocusing=min(length(i_refocusing),1);
1571 ktraj=zeros(size(gw));
1572 k=[0;0;0];
1573 for i=1:(length(i_periods)-1)
1574 %k=k+gw(:,i)*obj.gradRasterTime;
1575 i_period_end=(i_periods(i+1)-1);
1576 % here we use a trick to add current k value to the cumsum()
1577 k_period=cumsum([k,gw(:,i_periods(i):i_period_end)*obj.gradRasterTime],2);
1578 ktraj(:,i_periods(i):i_period_end)=k_period(:,2:end); % remove the first 'dummy' sample (see the trick above)
1579 k=k_period(:,end);
1580 if ii_next_excitation>0 && i_excitation(ii_next_excitation)==i_period_end
1581 k(:)=0;
1582 ktraj(:,i_period_end)=NaN; % we use NaN-s to mark the excitation point, they interrupt the plots
1583 ii_next_excitation = min(length(i_excitation),ii_next_excitation+1);
1584 end
1585 if ii_next_refocusing>0 && i_refocusing(ii_next_refocusing)==i_period_end
1586 k=-k;
1587 ii_next_refocusing = min(length(i_refocusing),ii_next_refocusing+1);
1588 end
1589 end
1591 % now calculate the k-space positions at the ADC time points
1592 % sample the k-space positions at the ADC time points
1593 ktraj_adc=interp1((1:(size(ktraj,2)))*obj.gradRasterTime, ktraj', ktime)';
1594 t_adc=ktime; % we now also return the sampling time points
1595 end
1597 function labels = evalLabels(obj, varargin)
1598 %Evaluate Label values of the entire sequence or its part
1599 % evalLabels(seqObj) Returns the label values at the end of
1600 % the sequence. Return value of the function is the structure
1601 % 'labels' with fields named after the labels used in the
1602 % sequence. Only the fiels corresponding to the lables
1603 % actually used are created.
1604 %
1605 % evalLabels(...,'blockRange',[first last]) Evaluate label
1606 % values starting from the first specified block to the last
1607 % one.
1608 %
1609 % evalLabels(...,'init',labels_struct) Evaluate labels
1610 % assuming the initial values from 'labels_struct'. Useful if
1611 % evaluating labes block-by-block.
1612 %
1613 % evalLabels(...,'evolution',<flag>) Evaluate labels and
1614 % return the evolution depending on the provided <flag>,
1615 % which is one of 'none','adc','label','blocks': 'blocks'
1616 % means return label values for all blocks; 'adc' - only for
1617 % blocks containig ADC objects; 'label' - only for blocks
1618 % where labels are manipulated.
1619 %
1620 validEvolutionValues = {'none','adc','label','blocks'};
1621 persistent parser
1622 if isempty(parser)
1623 parser = inputParser;
1624 parser.FunctionName = 'evalLabels';
1625 parser.addParamValue('blockRange',[1 inf],@(x)(isnumeric(x) && length(x)==2));
1626 parser.addParamValue('init',struct([]),@(x)(isempty(x) || isstruct(x)));
1627 parser.addParamValue('evolution','none',@(x)any(validatestring(x,validEvolutionValues)));
1628 end
1629 parse(parser,varargin{:});
1630 opt = parser.Results;
1632 if isempty(opt.init)
1633 labels=struct();
1634 else
1635 labels=opt.init;
1636 end
1638 %if ~strcmp(opt.evolution,'none')
1639 label_evol={};
1640 %end
1642 if ~isfinite(opt.blockRange(2))
1643 opt.blockRange(2)=length(obj.blockEvents);
1644 end
1646 for iB=opt.blockRange(1):opt.blockRange(2)
1647 block = obj.getBlock(iB);
1648 if isfield(block,'label') %current block has labels
1649 for i=1:length(block.label)
1650 if strcmp(block.label(i).type,'labelinc')
1651 if ~isfield(labels,block.label(i).label)
1652 labels.(block.label(i).label)=0;
1653 end
1654 labels.(block.label(i).label)=...
1655 labels.(block.label(i).label)+block.label(i).value;
1656 else
1657 labels.(block.label(i).label)=block.label(i).value;
1658 end
1659 end
1660 if strcmp(opt.evolution,'label')
1661 label_evol{end+1}=labels;
1662 end
1663 end
1664 if strcmp(opt.evolution,'blocks') || ...
1665 (strcmp(opt.evolution,'adc') && ~isempty(block.adc))
1666 label_evol{end+1}=labels;
1667 end
1668 end
1669 n=length(label_evol);
1670 if n>1 && ~strcmp(opt.evolution,'none')
1671 % convert cell array of structures to a structure of arrays
1672 % %l = cell2mat(label_evol);
1673 % l = [label_evol{:}]; % step1: convert to array of structures
1674 % f = fields(label_store);
1675 % a = cell(2,length(f)); % step2: prepare argumet
1676 % for i=1:length(f)
1677 % a{1,i}=f{i};
1678 % a{2,i}=[l.(f{i})];
1679 % end
1680 % label_store=struct(a{:}); % step3: create the final structure
1681 f = fieldnames(labels);
1682 for i=1:length(f)
1683 labels.(f{i})=zeros(1,n);
1684 for j=1:n
1685 if isfield(label_evol{j},f{i})
1686 labels.(f{i})(j)=label_evol{j}.(f{i});
1687 end
1688 end
1689 end
1690 end
1691 end
1693 function sp = plot(obj, varargin)
1694 %plot Plot the sequence in a new figure.
1695 % plot(seqObj) Plot the sequence
1696 %
1697 % Generates "classical" Pulseq 6-panel sequence plot. The
1698 % panels are "ADC", "RF magnitude", "RF phase", and three
1699 % gradient axes. ADC panel additionally shows digital output
1700 % pulses (a.k.a. triggers) as green diamonds with a line
1701 % indicatin the duration of the trigger pulse and input
1702 % triggers (e.g. cardiac) as triangles with a dot at the
1703 % center. ADC panal optionally also shows evolution of data
1704 % labels. RF phase panel shows RF phase at the center of the
1705 % RF pulse as a cross.
1706 %
1707 % plot(...,'timeRange',[start stop]) Plot the sequence
1708 % between the times specified by start and stop.
1709 %
1710 % plot(...,'blockRange',[first last]) Plot the sequence
1711 % starting from the first specified block to the last one.
1712 %
1713 % plot(...,'timeDisp',unit) Display time in:
1714 % 's', 'ms' or 'us'.
1715 %
1716 % plot(...,'label','LIN,REP') Plot label values for ADC events:
1717 % in this example for LIN and REP labels; other valid labes are
1718 % accepted as a comma-separated list.
1719 %
1720 % plot(...,'showBlocks',1) Plot grid and tick labels at the
1721 % block boundaries. Accepts a numeric or a boolean parameter.
1722 %
1723 % plot(...,'stacked',1) Rearrange the plots such they are vertically
1724 % stacked and share the same x-axis. Accepts a numeric or a boolean
1725 % parameter.
1726 %
1727 % plot(...,'showGuides',1) How dynamic hairline guides that follow
1728 % the data cursor to help verifying event alignment. Accepts a
1729 % numeric or a boolean parameter.
1730 %
1731 % f=plot(...) Return the new figure handle.
1732 %
1734 if nargout == 1
1735 sp = mr.aux.SeqPlot(obj, varargin{:});
1736 else
1737 mr.aux.SeqPlot(obj, varargin{:});
1738 end
1740 end
1742 function sp = paperPlot(obj, varargin)
1743 %paperPlot Plot the sequence in a stzle similar to that used in
1744 % scientific papers.
1745 % paperPlot(seqObj) Plot the sequence
1746 %
1747 % paperPlot(...,'blockRange',[first last]) Plot the sequence
1748 % starting from the first specified block to the last one.
1749 %
1750 % paperPlot(...,'lineWidth', w) Plot the sequence
1751 % using the specified line width.
1752 %
1753 % paperPlot(...,'axesColor', w) Plot the sequence
1754 % using the specified color for the horisontal axes.
1755 %
1756 % paperPlot(...,'rfColor', w) Plot the sequence
1757 % using the specified color for the RF and ADC events.
1758 %
1759 % paperPlot(...,'gxColor', w) Plot the sequence
1760 % using the specified color for the X gradients.
1761 %
1762 % paperPlot(...,'gyColor', w) Plot the sequence
1763 % using the specified color for the Y gradients.
1764 %
1765 % paperPlot(...,'gzColor', w) Plot the sequence
1766 % using the specified color for the Z gradients.
1767 %
1768 % paperPlot(...,'rfPlot', <'abs', 'real', 'imag'>) Plot the
1769 % RF pulses as the magnitude or real or imaginary part.
1770 %
1771 % Color parameters can be provided as a common color names,
1772 % e.g. 'red', 'blue', 'black', character strings starting
1773 % from '#' followed by a hexadecimal RGB values ranging from
1774 % 00 to ff or as an 1x3 vector of doubles ranging from 0 to 1
1775 % containing RGB values.
1776 %
1777 % f=paperPlot(...) Return the new figure handle.
1778 %
1780 function c=my_validatecolor(c)
1781 try
1782 c=validatecolor(c);
1783 catch
1784 c=[];
1785 end
1786 end
1787 validRfPlotValues = {'abs','real','imag'};
1788 persistent parser
1789 if isempty(parser)
1790 parser = inputParser;
1791 parser.FunctionName = 'paperPlot';
1792 parser.addParamValue('blockRange',[1 inf],@(x)(isnumeric(x) && length(x)==2));
1793 parser.addParamValue('lineWidth',1.2,@(x)(isnumeric(x)));
1794 parser.addParamValue('axesColor',[0.5 0.5 0.5],@(x)~isempty(validatecolor(x)));
1795 parser.addParamValue('rfColor','black',@(x)~isempty(my_validatecolor(x)));
1796 parser.addParamValue('gxColor','blue',@(x)~isempty(my_validatecolor(x)));
1797 parser.addParamValue('gyColor','red',@(x)~isempty(my_validatecolor(x)));
1798 parser.addParamValue('gzColor',[0 0.5 0.3],@(x)~isempty(my_validatecolor(x)));
1799 parser.addParamValue('rfPlot','abs',@(x)any(validatestring(x,validPlotRfValues)));
1800 end
1801 parse(parser,varargin{:});
1802 opt = parser.Results;
1804 if mr.aux.isOctave()
1805 warning('Function paperPlot() does not (yet) work on Octave.');
1806 return;
1807 end
1809 lw=opt.lineWidth;
1810 axes_clr=opt.axesColor;
1812 blockRange=opt.blockRange;
1813 if ~isfinite(blockRange(2))
1814 blockRange(2)=length(obj.blockDurations);
1815 end
1817 [wave_data,~,~,t_adc]=obj.waveforms_and_times(true,blockRange); % also export RF
1819 gwm=max(abs([wave_data{1:3}]'));
1820 rfm=max(abs([wave_data{4}]'));
1821 gwm(1)=max(gwm(1),t_adc(end));
1823 % remove horizontal lines with 0s. we detect 0 0 and insert a NaN in between
1824 for i=1:4
1825 j=size(wave_data{i},2);
1826 % this worked but was very slow...
1827 %while j>1
1828 % if wave_data{i}(2,j)==0 && wave_data{i}(2,j-1)==0
1829 % wave_data{i}(:,j:end+1)=[ [0.5*(wave_data{i}(1,j-1)+wave_data{i}(1,j));NaN] wave_data{i}(:,j:end)];
1830 % end
1831 % j=j-1;
1832 %end
1833 iInserts=find((wave_data{i}(2,1:end-1)==0) .* (wave_data{i}(2,2:end)==0));
1834 if ~isempty(iInserts)
1835 newWave=zeros(2,size(wave_data{i},2)+length(iInserts));
1836 c=1;
1837 for j=1:length(iInserts)
1838 newWave(:,(c+j-1):(iInserts(j)+j-1))=wave_data{i}(:,c:iInserts(j));
1839 newWave(:,iInserts(j)+j)=[0.5*(wave_data{i}(1,iInserts(j))+wave_data{i}(1,iInserts(j)+1));NaN];
1840 c=iInserts(j)+1;
1841 end
1842 newWave(:,(c+length(iInserts)):end)=wave_data{i}(:,c:end);
1843 wave_data{i}=newWave;
1844 end
1845 end
1847 f=figure;
1848 %f=colordef(f,'white'); %Set color scheme
1850 f.Color='w'; %Set background color of figure window
1852 t = tiledlayout(4,1,'TileSpacing','none');
1853 ax=[];
1855 nexttile
1856 % plot the 'axis'
1857 plot([-0.01*gwm(1),1.01*gwm(1)],[0 0],'Color',axes_clr,'LineWidth',lw/5); hold on;
1858 % plot the RF waveform
1859 %wave_data{4}(2,wave_data{4}(2,:)==0)=NaN; % hide 0s
1860 switch opt.rfPlot
1861 case 'real'
1862 plot(wave_data{4}(1,:), real(wave_data{4}(2,:)),'Color',opt.rfColor,'LineWidth',lw);
1863 case 'imag'
1864 plot(wave_data{4}(1,:), imag(wave_data{4}(2,:)),'Color',opt.rfColor,'LineWidth',lw);
1865 otherwise
1866 plot(wave_data{4}(1,:), abs(wave_data{4}(2,:)),'Color',opt.rfColor,'LineWidth',lw);
1867 end
1869 % plot ADCs
1870 t_adc_x3=repmat(t_adc,[3 1]);
1871 y_adc_x3=repmat([0; rfm(2)/5; NaN],[1 length(t_adc)]);
1872 plot(t_adc_x3(:),y_adc_x3(:),'Color',opt.rfColor,'LineWidth',lw/4);
1874 xlim([-0.03*gwm(1),1.03*gwm(1)]);
1875 ylim([-1.03*rfm(2),1.03*rfm(2)]);
1876 set(gca, 'box','off','XTickLabel',[],'XTick',[],'YTickLabel',[],'YTick',[]);
1877 set(get(gca, 'XAxis'), 'Visible', 'off');
1878 set(get(gca, 'YAxis'), 'Visible', 'off');
1879 ax(end+1)=gca;
1881 nexttile
1882 % plot the 'axis'
1883 plot([-0.01*gwm(1),1.01*gwm(1)],[0 0],'Color',axes_clr,'LineWidth',lw/5); hold on;
1884 % plot the entire gradient waveforms
1885 plot(wave_data{3}(1,:), wave_data{3}(2,:),'Color',opt.gzColor,'LineWidth',lw);
1887 xlim([-0.03*gwm(1),1.03*gwm(1)]);
1888 ylim([-1.03*gwm(2),1.03*gwm(2)]);
1889 set(gca, 'box','off','XTickLabel',[],'XTick',[],'YTickLabel',[],'YTick',[]);
1890 set(get(gca, 'XAxis'), 'Visible', 'off');
1891 set(get(gca, 'YAxis'), 'Visible', 'off');
1892 ax(end+1)=gca;
1894 nexttile
1895 % plot the 'axis'
1896 plot([-0.01*gwm(1),1.01*gwm(1)],[0 0],'Color',axes_clr,'LineWidth',lw/5); hold on;
1897 % plot the entire gradient waveforms
1898 plot(wave_data{2}(1,:), wave_data{2}(2,:),'Color',opt.gyColor,'LineWidth',lw);
1900 xlim([-0.03*gwm(1),1.03*gwm(1)]);
1901 ylim([-1.03*gwm(2),1.03*gwm(2)]);
1902 set(gca, 'box','off','XTickLabel',[],'XTick',[],'YTickLabel',[],'YTick',[]);
1903 set(get(gca, 'XAxis'), 'Visible', 'off');
1904 set(get(gca, 'YAxis'), 'Visible', 'off');
1905 ax(end+1)=gca;
1907 nexttile
1908 % plot the 'axis'
1909 plot([-0.01*gwm(1),1.01*gwm(1)],[0 0],'Color',axes_clr,'LineWidth',lw/5); hold on;
1910 % plot the entire gradient waveforms
1911 plot(wave_data{1}(1,:), wave_data{1}(2,:),'Color',opt.gxColor,'LineWidth',lw);
1913 xlim([-0.03*gwm(1),1.03*gwm(1)]);
1914 ylim([-1.03*gwm(2),1.03*gwm(2)]);
1915 set(gca, 'box','off','XTickLabel',[],'XTick',[],'YTickLabel',[],'YTick',[]);
1916 set(get(gca, 'XAxis'), 'Visible', 'off');
1917 set(get(gca, 'YAxis'), 'Visible', 'off');
1918 ax(end+1)=gca;
1920 % link zooming on the time axis
1921 linkaxes(ax(:),'x')
1923 if nargout == 1
1924 sp = f;
1925 end
1926 end
1928 function [wave_data, tfp_excitation, tfp_refocusing, t_adc, fp_adc, pm_adc]=waveforms_and_times(obj, appendRF, blockRange)
1929 % waveforms_and_times()
1930 % Decompress the entire gradient waveform
1931 % Returns gradient wave forms as a cell array with
1932 % gradient_axes (typically 3) dimension; each cell contains
1933 % time points and the correspndig gradient amplitude values.
1934 % Additional return values are time points of excitations,
1935 % refocusings and ADC sampling points.
1936 % If the optionl parameter 'appendRF' is set to true the RF
1937 % wave shapes are appended after the gradients
1938 % Optional output parameters: tfp_excitation contains time
1939 % moments, frequency and phase offsets of the excitation RF
1940 % pulses (similar for tfp_refocusing); t_adc contains times
1941 % of all ADC sample points; fp_adc contains frequency and
1942 % phase offsets of each ADC object (not sample); pm_adc
1943 % contains phase modulation of every adc sample beyond the
1944 % data stored in fp_adc (phaseModulation fields of v1.5.0).
1945 % TODO: return RF frequency offsets and RF waveforms and t_preparing (once its available)
1947 if nargin < 3
1948 blockRange=[1, length(obj.blockEvents)];
1949 else
1950 if length(blockRange)~=2
1951 error('parameter ''blockRange'' must contain exactly two numbers: first and last blocks from the range to consider');
1952 end
1953 end
1955 if nargin < 2
1956 appendRF=false;
1957 end
1959 grad_channels=3;
1960 gradChannels={'gx','gy','gz'}; % FIXME: this is not OK for matrix gradient systems
1962 t0=0;
1963 t0_n=0;
1965 numBlocks=blockRange(2)-blockRange(1)+1;
1967 % collect the shape pieces into a cell array
1968 if appendRF
1969 shape_channels=length(gradChannels)+1; % the last "channel" is RF
1970 else
1971 shape_channels=length(gradChannels);
1972 end
1973 shape_pieces=cell(shape_channels,numBlocks);
1974 % also collect RF and ADC timing data
1975 % t_excitation t_refocusing t_adc
1976 tfp_excitation=[];
1977 tfp_refocusing=[];
1978 t_adc=[];
1979 fp_adc=[];
1980 pm_adc=[];
1981 %block_durations=zeros(1,numBlocks);
1982 curr_dur=0;
1983 iP=0;
1984 out_len=zeros(1,shape_channels); % the last "channel" is RF
1985 for iBc=blockRange(1):blockRange(2)
1986 block = obj.getBlock(iBc);
1988 if isfield(block,'rotation')
1989 % apply the rotation to the current block and restore the block structure
1990 c=mr.rotate3D(block.rotation.rotQuaternion,block,'system',obj.sys);
1991 for i=1:3
1992 block.(gradChannels{i})=[];
1993 end
1994 for i=1:length(c)
1995 if isstruct(c{i}) && isfield(c{i},'type') && isfield(c{i},'channel')
1996 block.(['g' c{i}.channel])=c{i};
1997 end
1998 end
1999 end
2001 iP=iP+1;
2002 for j=1:length(gradChannels)
2003 grad=block.(gradChannels{j});
2004 if ~isempty(block.(gradChannels{j}))
2005 if strcmp(grad.type,'grad')
2006 % check if we have an extended trapezoid or an arbitrary gradient
2007 % on a regular raster. Arbitrary gradient on a pure centers raster
2008 % (shifted by 0.5) needs special processing
2009 tt_rast=grad.tt/obj.gradRasterTime;
2010 if all(abs(tt_rast-((1:length(tt_rast))-0.5)')<1e-6)
2011 % arbitrary gradient on a centers raster (no oversampling)
2012 % restore shape: if we had a trapezoid converted to shape we
2013 % have to find the "corners" and we can eliminate internal
2014 % samples on the straight segments but first we have to
2015 % restore samples on the edges of the gradient raster
2016 % intervals - for that we need the first sample
2018 [tt_chg, waveform_chg] = mr.restoreAdditionalShapeSamples(grad.tt,grad.waveform,grad.first,grad.last,obj.gradRasterTime,iBc);
2020 out_len(j)=out_len(j)+length(tt_chg);
2021 shape_pieces{j,iP}=[curr_dur+grad.delay+tt_chg; waveform_chg];%curr_dur+grad.delay+tgc;
2022 else
2023 % extended trapezoid or sampled gradient with oversampling (the easy case!)
2024 % the only caveat is that we need to add the first and last poins to the
2025 % shape in case of the oversampled rasterized gradient
2026 if abs(tt_rast(1)-0.5)<1e-6 % rasterized gradient's first sample is always on half-raster, extended trapezoid is always on a raster edje
2027 out_len(j)=out_len(j)+length(grad.tt)+2;
2028 shape_pieces{j,iP}=[curr_dur+grad.delay+[0 grad.tt' grad.shape_dur]; [grad.first grad.waveform' grad.last]];
2029 else
2030 out_len(j)=out_len(j)+length(grad.tt);
2031 shape_pieces{j,iP}=[curr_dur+grad.delay+grad.tt'; grad.waveform'];
2032 end
2033 end
2034 else
2035 if (abs(grad.flatTime)>eps) % interp1 gets confused by triangular gradients (repeating sample)
2036 out_len(j)=out_len(j)+4;
2037 shape_pieces{j,iP}=[
2038 curr_dur+grad.delay+cumsum([0 grad.riseTime grad.flatTime grad.fallTime]);...
2039 grad.amplitude*[0 1 1 0]];
2040 else
2041 if (abs(grad.riseTime)>eps && abs(grad.fallTime)>eps) % we skip 'empty' gradients
2042 out_len(j)=out_len(j)+3;
2043 shape_pieces{j,iP}=[
2044 curr_dur+grad.delay+cumsum([0 grad.riseTime grad.fallTime]);...
2045 grad.amplitude*[0 1 0]];
2046 else
2047 if abs(grad.amplitude)>eps
2048 warning('''empty'' gradient with non-zero magnitude detected in block %d',iBc);
2049 end
2050 end
2051 end
2052 end
2053 end
2054 end
2055 if ~isempty(block.rf)
2056 rf=block.rf;
2057 tc=mr.calcRfCenter(rf);
2058 t=rf.delay+tc;
2059 full_freqOffset=rf.freqOffset+rf.freqPPM*1e-6*obj.sys.gamma*obj.sys.B0;
2060 full_phaseOffset=rf.phaseOffset+rf.phasePPM*1e-6*obj.sys.gamma*obj.sys.B0;
2061 if (~isfield(rf,'use') || strcmp(rf.use,'excitation') || strcmp(rf.use,'undefined'))
2062 tfp_excitation(:,end+1) = [curr_dur+t; full_freqOffset; full_phaseOffset+2*pi*full_freqOffset*tc];
2063 elseif strcmp(rf.use,'refocusing')
2064 tfp_refocusing(:,end+1) = [curr_dur+t; full_freqOffset; full_phaseOffset+2*pi*full_freqOffset*tc];
2065 end
2066 if appendRF
2067 pre=[];
2068 post=[];
2069 if abs(rf.signal(1))>0
2070 pre=[curr_dur+rf.delay+rf.t(1)-eps;0];
2071 end
2072 if abs(rf.signal(end))>0
2073 post=[curr_dur+rf.delay+rf.t(end)+eps;0];
2074 end
2075% pre=[curr_dur+rf.delay+rf.t(1)-eps;NaN];
2076% post=[curr_dur+rf.delay+rf.t(end)+eps;NaN];
2077 out_len(end)=out_len(j)+length(rf.t)+size(pre,2)+size(post,2);
2078 shape_pieces{end,iP}=[pre [curr_dur+rf.delay+rf.t.'; (rf.signal.*exp(1i*(full_phaseOffset+2*pi*full_freqOffset*rf.t))).'] post];
2079 end
2080 end
2081 if ~isempty(block.adc)
2082 ta=block.adc.dwell*((0:(block.adc.numSamples-1))+0.5); % according to the information from Klaus Scheffler and indirectly from Siemens this is the present convention (the samples are shifted by 0.5 dwell) % according to the information from Klaus Scheffler and indirectly from Siemens this is the present convention (the samples are shifted by 0.5 dwell)
2083 n_adc_samples=length(t_adc);
2084 t_adc((end+1):(end+block.adc.numSamples)) = ta + block.adc.delay + curr_dur;
2085 full_freqOffset=block.adc.freqOffset+block.adc.freqPPM*1e-6*obj.sys.gamma*obj.sys.B0;
2086 full_phaseOffset=block.adc.phaseOffset+block.adc.phasePPM*1e-6*obj.sys.gamma*obj.sys.B0;
2087 if isempty(block.adc.phaseModulation)
2088 block.adc.phaseModulation=0;
2089 if nargout>=6
2090 pm_adc((n_adc_samples+1):(n_adc_samples+block.adc.numSamples))=zeros(1,block.adc.numSamples);
2091 end
2092 else
2093 if nargout>=6
2094 pm_adc((n_adc_samples+1):(n_adc_samples+block.adc.numSamples))=block.adc.phaseModulation;
2095 end
2096 end
2097 fp_adc(:,(end+1):(end+block.adc.numSamples)) = [full_freqOffset*ones(1,block.adc.numSamples); full_phaseOffset+block.adc.phaseModulation+full_freqOffset*ta];
2098 end
2099 curr_dur=curr_dur+obj.blockDurations(iBc);%mr.calcDuration(block);
2100 end
2102 % collect wave data
2103 wave_data=cell(1,shape_channels);
2104 for j=1:shape_channels
2105 wave_data{j}=zeros(2,out_len(j));
2106 end
2107 wave_cnt=zeros(1,shape_channels);
2108 curr_dur=0;
2109 for iP=1:numBlocks
2110 for j=1:shape_channels
2111 if ~isempty(shape_pieces{j,iP})
2112 wave_data_local=shape_pieces{j,iP};
2113 len=size(wave_data_local,2);
2114 if wave_cnt(j)~=0 && wave_data{j}(1,wave_cnt(j))+obj.gradRasterTime < wave_data_local(1,1)
2115 if wave_data{j}(2,wave_cnt(j))~=0
2116 if abs(wave_data{j}(2,wave_cnt(j)))>1e-6 % todo: real physical tolarance for gradient amplitudes
2117 warning('waveforms_and_times(): forcing ramp-down from a non-zero gradient sample on axis %d at t=%d us \ncheck your sequence, some calculations are possibly wrong. If using mr.makeArbitraryGrad() consider using explicit values for ''first'' and ''last'' and setting them correctly.', j, round(1e6*wave_data{j}(1,wave_cnt(j))));
2118 wave_data{j}(:,wave_cnt(j)+1)=[wave_data{j}(1,wave_cnt(j))+obj.gradRasterTime/2; 0]; % this is likely to cause memory reallocations
2119 wave_cnt(j)=wave_cnt(j)+1;
2120 else
2121 % we are within the tolorance, just set it to 0 quietly
2122 wave_data{j}(2,wave_cnt(j))=0.0;
2123 end
2124 end
2125 if wave_data_local(2,1)~=0
2126 if abs(wave_data_local(2,1))>1e-6 % todo: real physical tolarance for gradient amplitudes
2127 warning('waveforms_and_times(): forcing ramp-up to a non-zero gradient sample on axis %d at t=%d us \ncheck your sequence, some calculations are probably wrong. If using mr.makeArbitraryGrad() consider using explicit values for ''first'' and ''last'' and setting them correctly.', j, round(1e6*wave_data_local(1,1)));
2128 wave_data_local=[[wave_data_local(1,1)-obj.gradRasterTime/2; 0] wave_data_local]; % this is likely to cause memory reallocations also later on
2129 len=len+1;
2130 else
2131 % we are wihin the tolorance, just set it to 0 quietly
2132 wave_data_local(2,1)=0.0;
2133 end
2134 end
2135 end
2136 if wave_cnt(j)==0 || wave_data{j}(1,wave_cnt(j))<wave_data_local(1,1)
2137 wave_data{j}(:,wave_cnt(j)+(1:len))=wave_data_local;
2138 wave_cnt(j)=wave_cnt(j)+len;
2139 else
2140 if (wave_data_local(1,1)<wave_data{j}(1,wave_cnt(j))-1e-9) % TODO consistent time tolerance
2141 warning('Warning: looks like rounding errors for some elements exceed the acceptable tolerance!\n');
2142 end
2143 [~,d]=find(wave_data_local(1,:)>wave_data{j}(1,wave_cnt(j)),1);
2144 wave_data{j}(:,wave_cnt(j)+(1:(len-d+1)))=wave_data_local(:,d:end);
2145 wave_cnt(j)=wave_cnt(j)+len-d+1;
2146 end
2147 end
2148 end
2149 end
2150 for j=1:shape_channels
2151 if any(diff(wave_data{j}(1,1:wave_cnt(j)))<=0.0) %&& ... % quick pre-check whether the time vector is monotonously increasing to avoid too often unique() calls
2152 %wave_cnt(j)~=length(unique(wave_data{j}(1,1:wave_cnt(j))))
2153 warning('Warning: not all elements of the generated time vector are unique and sorted in accending order!');
2154 end
2155 end
2158 % trim the output data
2159 for j=1:shape_channels
2160 if wave_cnt(j)<size(wave_data{j},2)
2161 wave_data{j}(:,(wave_cnt(j)+1):end)=[];
2162 end
2163 end
2165% % convert wave data to piecewise polynomials
2166% wave_pp=cell(1,length(gradChannels));
2167% for j=1:length(gradChannels)
2168% if (wave_cnt(j)<=0)
2169% continue;
2170% end
2171% if ~all(isfinite(wave_data{j}(:)))
2172% fprintf('Warning: not all elements of the generated waveform are finite!\n');
2173% end
2174% wave_pp{j} = interp1(wave_data{j}(1,1:wave_cnt(j)),wave_data{j}(2,1:wave_cnt(j)),'linear','pp');
2175% end
2176 end
2178 function [ktraj_adc, t_adc, ktraj, t_ktraj, t_excitation, t_refocusing, slicepos, t_slicepos, gw_pp, pm_adc] = calculateKspacePP(obj, varargin)
2179 % calculate the k-space trajectory of the entire pulse sequence
2180 % using piecewise-polynomial gradient wave representation
2181 % which is much faster for simple shapes and large delays
2182 % optional parameter 'trajectory_delay' sets the compensation
2183 % factor to align ADC and gradients in the reconstruction
2184 % optional parameter 'gradient_offset' allows to simulate
2185 % background gradients or verifz spin-echo conditions
2186 % Return values: ktraj_adc, t_adc, ktraj, t_ktraj,
2187 % t_excitation, t_refocusing, slicepos
2189 persistent parser
2190 if isempty(parser)
2191 parser = inputParser;
2192 parser.FunctionName = 'calculateKspacePP';
2193 parser.addParamValue('trajectory_delay',0,@(x)(isnumeric(x)));
2194 parser.addParamValue('gradient_offset',0,@(x)(isnumeric(x)));
2195 parser.addParamValue('blockRange',[1 inf],@(x)(isnumeric(x) && length(x)==2));
2196 parser.addParamValue('externalWaveformsAndTimes',struct([]),@(x)(isstruct(x)));
2197 end
2198 parse(parser,varargin{:});
2199 opt = parser.Results;
2201 if any(abs(opt.trajectory_delay)>100e-6)
2202 warning('trajectory delay of (%s) us is suspiciously high',num2str(opt.trajectory_delay*1e6));
2203 end
2205 blockRange=opt.blockRange;
2206 if blockRange(1)<1
2207 blockRange(1)=1;
2208 end
2209 if ~isfinite(blockRange(2))
2210 blockRange(2)=length(obj.blockDurations);
2211 end
2213 total_duration=sum(obj.blockDurations(blockRange(1):blockRange(2)));
2215 if isempty(opt.externalWaveformsAndTimes)
2216 if nargout>=10
2217 [gw_data, tfp_excitation, tfp_refocusing, t_adc, ~,pm_adc]=obj.waveforms_and_times(false,blockRange);
2218 else
2219 [gw_data, tfp_excitation, tfp_refocusing, t_adc]=obj.waveforms_and_times(false,blockRange);
2220 end
2221 else
2222 gw_data=opt.externalWaveformsAndTimes.gw_data;
2223 tfp_excitation=opt.externalWaveformsAndTimes.tfp_excitation;
2224 tfp_refocusing=opt.externalWaveformsAndTimes.tfp_refocusing;
2225 t_adc=opt.externalWaveformsAndTimes.t_adc;
2226 if isfield(opt.externalWaveformsAndTimes, 'pm_adc')
2227 pm_adc=opt.externalWaveformsAndTimes.pm_adc;
2228 else
2229 pm_adc=[];
2230 end
2231 % how do we verify that the total_duration is correct???
2232 end
2234 ng=length(gw_data);
2235 % gradient delay handling
2236 if length(opt.trajectory_delay)==1
2237 gradient_delays(1:ng)=opt.trajectory_delay;
2238 else
2239 assert(length(opt.trajectory_delay)==ng); % we need to have the same number of gradient channels
2240 gradient_delays=opt.trajectory_delay;
2241 end
2242 % gradient offset handling
2243 if length(opt.gradient_offset)==1
2244 gradient_offset(1:ng)=opt.gradient_offset;
2245 else
2246 assert(length(opt.gradient_offset)==ng); % we need to have the same number of gradient channels
2247 gradient_offset=opt.gradient_offset;
2248 end
2250 % convert wave data to piecewise polynomials
2251 gw_pp=cell(1,ng);
2252 for j=1:ng
2253 wave_cnt=size(gw_data{j},2);
2254 if wave_cnt==0
2255 if abs(gradient_offset(j))<=eps % gradient offset support, part 1
2256 continue;
2257 else
2258 gw=[0,total_duration; 0, 0];
2259 end
2260 else
2261 gw=gw_data{j};
2262 end
2263 % now gw contains the wave form for the current axis
2264 if abs(gradient_delays(j))>eps
2265 gw(1,:)=gw(1,:)-gradient_delays(j); % (anisotropic) gradient delay support
2266 end
2267 if ~all(isfinite(gw(:)))
2268 fprintf('Warning: not all elements of the generated waveform are finite!\n');
2269 end
2270 teps=1e-12; % eps is too small and may go lost due to rounding errors (e.g. total_duration+eps==total_duration)
2271 if gw(1,1)>0 && gw(1,end) < total_duration
2272 gw=[ [-teps gw(1,1)-teps;0 0] gw [gw(1,end)+teps total_duration+teps;0 0] ]; % we need these "eps" terms to avoid integration errors over extended periods of time
2273 elseif gw(1,1)>0
2274 gw=[ [-teps gw(1,1)-teps;0 0] gw ]; % we need these "eps" terms to avoid integration errors over extended periods of time
2275 elseif gw(1,end) < total_duration
2276 gw=[ gw [gw(1,end)+teps total_duration+teps;0 0] ]; % we need these "eps" terms to avoid integration errors over extended periods of time
2277 end
2278 %
2279 if abs(gradient_offset(j))>eps
2280 gw(2,:)=gw(2,:)+gradient_offset(j); % gradient offset support, part 2
2281 end
2282 gw_pp{j} = interp1(gw(1,:),gw(2,:),'linear','pp');
2283 end
2285 % calculate slice positions. for now we entirely rely on the
2286 % excitation -- ignoring complicated interleaved refocused sequences
2287 if ~isempty(tfp_excitation)
2288 slicepos=zeros(length(gw_data),size(tfp_excitation,2)); % position in x,y,z;
2289 for j=1:length(gw_data)
2290 if isempty(gw_pp{j})
2291 slicepos(j,:)=NaN;
2292 else
2293 slicepos(j,:)=tfp_excitation(2,:)./ppval(gw_pp{j},tfp_excitation(1,:));
2294 end
2295 end
2296 slicepos(~isfinite(slicepos))=0; % reset undefined to 0 (or is NaN better?)
2297 t_slicepos=tfp_excitation(1,:);
2298 else
2299 slicepos=[];
2300 t_slicepos=[];
2301 end
2303 %t_adc = t_adc + opt.trajectory_delay;
2304 % this was wrong because it dod not shift RF events (which are
2305 % intrinsically well-synchronized) and did not allow for
2306 % anisotropic delays for different gradient axes. For the new
2307 % implementation see "gradient_delays" vector above
2309 % integrate waveforms as PPs to produce gadient moments
2310 gm_pp=cell(1,ng);
2311 tc = {};
2312 for i=1:ng
2313 if isempty(gw_pp{i})
2314 continue;
2315 end
2316 if mr.aux.isOctave()
2317 gm_pp{i}=ppint(gw_pp{i});
2318 else
2319 gm_pp{i}=fnint(gw_pp{i});
2320 end
2321 tc{end+1}=gm_pp{i}.breaks;
2322 % "sample" ramps for display purposes otherwise piecewise-linear diplay (plot) fails (looks stupid)
2323 ii=find(abs(gm_pp{i}.coefs(:,1))>eps);
2324 if ~isempty(ii)
2325 tca=cell(1,length(ii));
2326 for j=1:length(ii)
2327 tca{j}=(floor(gm_pp{i}.breaks(ii(j))/obj.gradRasterTime):ceil(gm_pp{i}.breaks(ii(j)+1)/obj.gradRasterTime))*obj.gradRasterTime;
2328 end
2329 tc{end+1}=[tca{:}];
2330 end
2331 end
2332 %t = unique([tc{:}, 0, t_excitation-obj.gradRasterTime, t_excitation, t_refocusing, t_adc]);
2333 % we round to 100ns, otherwise unique() fails...
2335 if isempty(tfp_excitation)
2336 t_excitation=[];
2337 else
2338 t_excitation=tfp_excitation(1,:);
2339 end
2340 if isempty(tfp_refocusing)
2341 t_refocusing=[];
2342 else
2343 t_refocusing=tfp_refocusing(1,:);
2344 end
2346 tacc=1e-10; % temporal accuracy
2347 taccinv=1/tacc;
2348 t_ktraj = tacc*unique(round(taccinv*[tc{:}, 0, t_excitation-2*obj.rfRasterTime, t_excitation-obj.rfRasterTime, t_excitation, t_refocusing-obj.rfRasterTime, t_refocusing, t_adc, total_duration]));
2349 % % the "proper" matlab's function ismember() is slow and returns a bool array, but builtin('_ismemberhelper'...) is not compatible accross versions (known not to work on the windows 2021a version)
2350 %[~,i_excitation]=builtin('_ismemberhelper',tacc*round(taccinv*t_excitation),t_ktraj);
2351 %[~,i_refocusing]=builtin('_ismemberhelper',tacc*round(taccinv*t_refocusing),t_ktraj);
2352 %[~,i_adc]=builtin('_ismemberhelper',tacc*round(taccinv*t_adc),t_ktraj);
2353 % this is nother undocumented solution, see https://undocumentedmatlab.com/blog_old/ismembc-undocumented-helper-function
2354 if mr.aux.isOctave()
2355 [~,i_excitation]=ismember(tacc*round(taccinv*t_excitation),t_ktraj);
2356 [~,i_refocusing]=ismember(tacc*round(taccinv*t_refocusing),t_ktraj);
2357 [~,i_adc]=ismember(tacc*round(taccinv*t_adc),t_ktraj);
2358 else
2359 i_excitation=ismembc2(tacc*round(taccinv*t_excitation),t_ktraj);
2360 i_refocusing=ismembc2(tacc*round(taccinv*t_refocusing),t_ktraj);
2361 i_adc=ismembc2(tacc*round(taccinv*t_adc),t_ktraj);
2362 end
2363 %
2364 i_periods=unique([1, i_excitation, i_refocusing, length(t_ktraj)]);
2365 if ~isempty(i_excitation)
2366 ii_next_excitation=1;
2367 else
2368 ii_next_excitation=0;
2369 end
2370 if ~isempty(i_refocusing)
2371 ii_next_refocusing=1;
2372 else
2373 ii_next_refocusing=0;
2374 end
2375 ktraj=zeros(3, length(t_ktraj));
2376 for i=1:ng
2377 if isempty(gw_pp{i})
2378 continue;
2379 end
2380 %[~,it]=builtin('_ismemberhelper',[gm_pp{i}.breaks(1),gm_pp{i}.breaks(end)],t);
2381 %ktraj(i,it(1):it(2))=ppval(gm_pp{i},t(it(1):it(2)));
2382 it=find(t_ktraj>=tacc*round(taccinv*gm_pp{i}.breaks(1)) & t_ktraj<=tacc*round(taccinv*gm_pp{i}.breaks(end)));
2383 ktraj(i,it)=ppval(gm_pp{i},t_ktraj(it));
2384 if t_ktraj(it(end))<t_ktraj(end)
2385 ktraj(i,(it(end)+1):end)=ktraj(i,it(end));
2386 end
2387 end
2388 % convert gradient moments to k-space
2389 dk=-ktraj(:,1);%[0;0;0];
2390 for i=1:(length(i_periods)-1)
2391 i_period=i_periods(i);
2392 i_period_end=i_periods(i+1);
2393 if ii_next_excitation>0 && i_excitation(ii_next_excitation)==i_period
2394 if abs(t_ktraj(i_period)-t_excitation(ii_next_excitation))>tacc
2395 warning('abs(t_ktraj(i_period)-t_excitation(ii_next_excitation))<%g failed for ii_next_excitation=%d error=%g', tacc, ii_next_excitation, t_ktraj(i_period)-t_excitation(ii_next_excitation));
2396 end
2397 dk=-ktraj(:,i_period);
2398 if (i_period>1)
2399 ktraj(:,i_period-1)=NaN; % we use NaN-s to mark the excitation point, they interrupt the plots
2400 end
2401 ii_next_excitation = min(length(i_excitation),ii_next_excitation+1);
2402 elseif ii_next_refocusing>0 && i_refocusing(ii_next_refocusing)==i_period
2403 %dk=-ktraj(:,i_period);
2404 dk=-2*ktraj(:,i_period)-dk;
2405 %if (i_period>1)
2406 % ktraj(:,i_period-1)=NaN; % we use NaN-s to mark the excitation point, they interrupt the plots
2407 %end
2408 ii_next_refocusing = min(length(i_refocusing),ii_next_refocusing+1);
2409 end
2411 ktraj(:,i_period:(i_period_end-1))=ktraj(:,i_period:(i_period_end-1))+dk;
2412 end
2413 ktraj(:,i_period_end)=ktraj(:,i_period_end)+dk;
2414 ktraj_adc=ktraj(:,i_adc);
2415% % add first and last slicepos
2416% if t_slicepos(1)>0
2417% t_slicepos=[0 t_slicepos];
2418% slicepos=[ zeros(length(gw_data),1) slicepos ];
2419% end
2420% if t_slicepos(end)<t_ktraj(end)
2421% t_slicepos=[t_slicepos t_ktraj(end)];
2422% slicepos=[ slicepos slicepos(:,end)];
2423% end
2424 end
2426 function [mean_pwr, peak_pwr, rf_rms, total_energy]=calcRfPower(obj, varargin)
2427 %calcRfPower : Calculate the relative** power of the RF pulse
2428 % Returns the (relative) energy of the pulse expressed in the units of
2429 % RF amplitude squared multiplied by time, e.g. in Pulseq these are
2430 % Hz * Hz * s = Hz. Sounds strange, but is true. Return
2431 % parameter mean_pwr is closely related to the relative** SAR.
2432 % Also returns peak power [Hz^2] and RMS B1 amplitude [Hz].
2433 % Optional parameter 'blockRange' can be used to specify the
2434 % part of the sequence for which the energy is calculated.
2435 % Optional parameter 'windowDuration' can be used to specify
2436 % the time window for the total_energy, mean_pwr and rf_rms
2437 % calculation. The values returned in this case are the
2438 % maximum values over all time windows. The time window is
2439 % rounded up to a certain number of complete blocks.
2440 % ** Note: the power and rf amplitude calculated by this function is
2441 % relative as it is calculated in units of Hz^2 or Hz. The rf amplitude
2442 % can be converted to T by dividing the resulting value by gamma.
2443 % Correspondingly, The power can be converted to mT^2*s by dividing
2444 % the given value by gamma^2. Nonetheless, the absolute SAR is related to
2445 % the electric field, so the further scaling coeficient is both tx-coil-
2446 % dependent (e.g. depends on the coil design) and also subject-dependent
2447 % (e.g. depends on the reference voltage).
2449 persistent parser
2450 if isempty(parser)
2451 parser = inputParser;
2452 parser.FunctionName = 'calcRfPower';
2453 parser.addParamValue('blockRange',[1 inf],@(x)(isnumeric(x) && length(x)==2));
2454 parser.addParamValue('windowDuration',NaN,@(x)(isnumeric(x)));
2455 end
2456 parse(parser,varargin{:});
2457 opt = parser.Results;
2459 if ~isfinite(opt.blockRange(2))
2460 opt.blockRange(2)=length(obj.blockEvents);
2461 end
2463 dur=0;
2464 total_energy=0;
2465 peak_pwr=0;
2466 rf_ms=0;
2468 windowOn= isfinite(opt.windowDuration);
2469 if windowOn
2470 blockBookkeeping = zeros(2,opt.blockRange(2)-opt.blockRange(1)+1);
2471 currentWindowDur=0.0;
2472 currentWindowStartIdx=opt.blockRange(1);
2473 total_energy_max=0.0;
2474 rf_ms_max=0.0;
2475 end
2477 for iBc=opt.blockRange(1):opt.blockRange(2)
2478 block = obj.getBlock(iBc);
2479 dur=dur+obj.blockDurations(iBc);
2481 if ~isempty(block.rf)
2482 rf=block.rf;
2483 [e,pp,rms]=mr.calcRfPower(rf);
2485 total_energy=total_energy+e;
2486 rf_ms=rf_ms+rms^2*rf.shape_dur;
2487 peak_pwr=max(peak_pwr,pp);
2489 if windowOn
2490 blockBookkeeping(:,iBc-opt.blockRange(1)+1)=[e,rms^2*rf.shape_dur];
2491 total_energy_max=max(total_energy_max,total_energy);
2492 rf_ms_max=max(rf_ms_max,rf_ms);
2493 end
2494 end
2495 % keep track of the window width and make it shorter if needed
2496 if windowOn
2497 currentWindowDur=currentWindowDur+obj.blockDurations(iBc);
2498 while currentWindowDur>opt.windowDuration
2499 % remove the front side of the window from total_energy and rf_ms
2500 total_energy=total_energy-blockBookkeeping(1,currentWindowStartIdx);
2501 rf_ms=rf_ms-blockBookkeeping(2,currentWindowStartIdx);
2502 currentWindowDur=currentWindowDur-obj.blockDurations(currentWindowStartIdx);
2503 currentWindowStartIdx=currentWindowStartIdx+1;
2504 end
2505 end
2506 end
2508 if windowOn
2509 total_energy=total_energy_max;
2510 mean_pwr=total_energy/opt.windowDuration; % we divide here by the nominal window duration, which may lead to a slight overestimation
2511 rf_rms=sqrt(rf_ms_max/opt.windowDuration); % we divide here by the nominal window duration, which may lead to a slight overestimation
2512 else
2513 mean_pwr=total_energy/dur;
2514 rf_rms=sqrt(rf_ms/dur);
2515 end
2516 end
2518 function applySoftDelay(obj, varargin)
2519 % applies soft delays to the sequence by modifying the block
2520 % durations of the respective blocks. Input parameters are
2521 % pairs of soft delays and values, whereas the soft delay is
2522 % identified by its string hint and the value is the duration
2523 % in seconds. Not all soft delays defined in the sequence need
2524 % to be specified. Examples:
2525 %
2526 % seq.applySoftDelay('TE',40e-3); % set TE to 40ms
2527 % seq.applySoftDelay('TE',50e-3,'TR',2); % set TE to 50ms and TR to 2 s
2528 %
2529 % See also: mr.makeSoftDelay()
2531 % parse arguments manually
2532 sdm_input=containers.Map('KeyType', 'char', 'ValueType', 'double');
2533 for i=1:2:length(varargin)
2534 if ~ischar(varargin{i})
2535 error('Argument at the position %d must be a character string ID of the soft delay',i);
2536 end
2537 if i>length(varargin) || ~isnumeric(varargin{i+1})
2538 error('Argument at the position %d must be the value of the soft delay ''%s''',i+1,varargin{i});
2539 end
2540 sdm_input(varargin{i})=varargin{i+1};
2541 end
2542 % go through all the blocks and update durations, at the same time
2543 % checking the consistency of the soft delays
2544 sdm_Str2numIDs=containers.Map('KeyType', 'char', 'ValueType', 'double');
2545 sdm_num2hint=containers.Map('KeyType', 'double', 'ValueType', 'char');
2546 sdm_warnings=containers.Map('KeyType', 'double', 'ValueType', 'logical');
2547 for iBc=1:length(obj.blockDurations)
2548 b = obj.getBlock(iBc);
2549 if isfield(b, 'softDelay') && ~isempty(b.softDelay)
2550 % check the numeric ID consistency
2551 if ~sdm_Str2numIDs.isKey(b.softDelay.hint)
2552 sdm_Str2numIDs(b.softDelay.hint)=b.softDelay.num;
2553 else
2554 if sdm_Str2numIDs(b.softDelay.hint)~=b.softDelay.num
2555 error('Soft delay in block %d with numeric ID %d and string hint ''%s'' is inconsistent with the previous occurences of the same string hint', iBc, b.softDelay.num, b.softDelay.hint);
2556 end
2557 end
2558 if ~sdm_num2hint.isKey(b.softDelay.num)
2559 sdm_num2hint(b.softDelay.num)=b.softDelay.hint;
2560 else
2561 if ~strcmp(sdm_num2hint(b.softDelay.num),b.softDelay.hint)
2562 error('Soft delay in block %d with numeric ID %d and string hint ''%s'' is inconsistent with the previous occurences of the same numeric ID', iBc, b.softDelay.num, b.softDelay.hint);
2563 end
2564 end
2565 if sdm_input.isKey(b.softDelay.hint)
2566 % calculate the new block duration
2567 new_dur_ru=(sdm_input(b.softDelay.hint)/b.softDelay.factor + b.softDelay.offset)/obj.sys.blockDurationRaster;
2568 new_dur=round(new_dur_ru)*obj.sys.blockDurationRaster;
2569 if abs(new_dur-new_dur_ru*obj.sys.blockDurationRaster)>0.5e-6 && ~sdm_warnings.isKey(b.softDelay.num)
2570 warning('Block duration for block %d, soft delay ''%s'', had to be substantially rounded to become aligned to the raster time. This warning is only displayed for the first block where it occurs.', iBc, b.softDelay.hint);
2571 sdm_warnings(b.softDelay.num)=true;
2572 end
2573 if new_dur<0
2574 error('Calculated new duration of the block %i, soft delay %s/%d is negative (%g s)', iBc, b.softDelay.hint, b.softDelay.num, new_dur);
2575 end
2576 obj.blockDurations(iBc)=new_dur;
2577 end
2578 end
2579 end
2580 % now check if there are some input soft delays which haven't been found in the sequence
2581 all_input_hints=sdm_input.keys;
2582 for i=1:length(all_input_hints)
2583 if ~sdm_Str2numIDs.isKey(all_input_hints{i})
2584 error('Specified soft delay ''%s'' does not exist in the sequence', all_input_hints{i});
2585 end
2586 end
2587 end
2589 function [easyStruct, errorReport, softDelayState] =getDefaultSoftDelayValues(obj)
2590 % go through all the blocks checking the consistency of the soft delays
2591 % the code below is copied from checkTiming; we should merge
2592 % the functionality eventually... TODO/FIXME
2593 errorReport={};
2594 softDelayState={};
2595 for iB=1:length(obj.blockDurations)
2596 b = obj.getBlock(iB);
2597 % check soft delays
2598 if isfield(b, 'softDelay') && ~isempty(b.softDelay)
2599 if b.softDelay.factor==0
2600 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' soft delay ' b.softDelay.hint '/' num2str(b.softDelay.num) ' has factor parameter of 0 which is invalid\n' ] };
2601 is_ok=false;
2602 end
2603 % calculate the default delay value based on the current block duration
2604 def_del=(obj.blockDurations(iB)-b.softDelay.offset)*b.softDelay.factor;
2605 if (b.softDelay.num>=0)
2606 % remember or check for consistency
2607 if length(softDelayState)<b.softDelay.num+1 || isempty(softDelayState{b.softDelay.num+1})
2608 softDelayState{b.softDelay.num+1}=struct('def',def_del,'hint',b.softDelay.hint, 'blk', iB, 'min', 0.0, 'max', +Inf);
2609 else
2610 if abs(def_del-softDelayState{b.softDelay.num+1}.def)>1e-7 % what is the reasonable threshold?
2611 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' soft delay ' b.softDelay.hint '/' num2str(b.softDelay.num) ': default duration derived from this block (' num2str(def_del*1e6) 'us) is inconsistent with the previous default (' num2str(softDelayState{b.softDelay.num+1}.def*1e6) 'us) that was derived from block ' num2str(softDelayState{b.softDelay.num+1}.blk) '\n' ] };
2612 is_ok=false;
2613 end
2614 if ~strcmp(b.softDelay.hint, softDelayState{b.softDelay.num+1}.hint)
2615 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' soft delay ' b.softDelay.hint '/' num2str(b.softDelay.num) ': soft delays with the same numeric ID are expected to share the same text hint but previous hint recorded in block ' num2str(softDelayState{b.softDelay.num+1}.blk) ' is ' softDelayState{b.softDelay.num+1}.hint '\n' ] };
2616 is_ok=false;
2617 end
2618 end
2619 % calculate the delay value that would make the block duration of 0, which correponds to min/max
2620 lim_del=(-b.softDelay.offset)*b.softDelay.factor;
2621 if b.softDelay.factor>0
2622 % lim_del corresponds to a minimum
2623 if lim_del>softDelayState{b.softDelay.num+1}.min
2624 softDelayState{b.softDelay.num+1}.min=lim_del;
2625 end
2626 else
2627 % lim_del corresponds to a maximum
2628 if lim_del<softDelayState{b.softDelay.num+1}.max
2629 softDelayState{b.softDelay.num+1}.max=lim_del;
2630 end
2631 end
2632 else
2633 errorReport = { errorReport{:}, [ ' Block:' num2str(iB) ' contains a soft delay ' b.softDelay.hint ' with an invalid numeric ID' num2str(b.softDelay.num) '\n' ] };
2634 is_ok=false;
2635 end
2636 end
2637 end
2638 % re-package softDelayState into easyStruct
2639 easyStruct=struct;
2640 for i=1:length(softDelayState)
2641 if isempty(softDelayState{i})
2642 warning('SoftDelay numeric ID %d is unused, we expect contiguous numbering of soft delays',i-1);
2643 continue;
2644 end
2645 if isfield(easyStruct,softDelayState{i}.hint)
2646 error('SoftDelay with numeric ID %d uses the rame hint ''%s'' as some previous SoftDelay',i-1,softDelayState{i}.hint);
2647 continue;
2648 end
2649 easyStruct.(softDelayState{i}.hint)=softDelayState{i}.def;
2650 end
2651 end
2653 function soundData=sound(obj, varargin)
2654 %sound()
2655 % "play out" the sequence through the system speaker and
2656 % return the sound data (if needed). Optional parameters are
2657 % 'blockRange', 'channelWeights', 'sampleRate',and
2658 % 'onlyProduceSoundData'. The latter skips the "playing out"
2659 % part. Sound data is 2xN array sampled at the provided
2660 % 'sampleRate' (default is CD quality 4411 Hz). The sound
2661 % vector is produced from the gradient waveforms with X and Y
2662 % mapped to challens 1 and 2, respectively, and Z split
2663 % between channels 1 ans 2. The output is then
2664 % Gauss-filtered to reduce high-frequency ringing and
2665 % normalized to 0.95.
2666 %
2668 persistent parser
2669 if isempty(parser)
2670 parser = inputParser;
2671 parser.FunctionName = 'sound';
2672 parser.addParamValue('blockRange',[1 inf],@(x)(isnumeric(x) && length(x)==2));
2673 parser.addParamValue('channelWeights',[1 1 1],@(x)(isnumeric(x) && length(x)==3));
2674 parser.addParamValue('onlyProduceSoundData',false,@(x)(islogical(x)));
2675 parser.addParamValue('sampleRate',44100,@(x)(isnumeric(x) && isscalar(x) && x>0));
2676 end
2677 parse(parser,varargin{:});
2678 opt = parser.Results;
2680 if ~isfinite(opt.blockRange(2))
2681 opt.blockRange(2)=length(obj.blockEvents);
2682 end
2684 gw_data=obj.waveforms_and_times(false,opt.blockRange);
2685 total_duration=sum(obj.blockDurations);
2687 sample_rate=opt.sampleRate; % default is 44100 Hz (CD quality)
2688 dwell_time=1/sample_rate;
2689 sound_length=floor(total_duration/dwell_time)+1;
2691 soundData(2,sound_length)=0; %preallocate
2693 if ~isempty(gw_data{1})
2694 soundData(1,:)=interp1(gw_data{1}(1,:),gw_data{1}(2,:)*opt.channelWeights(1),(0:(sound_length-1))*dwell_time,'linear',0);
2695 end
2696 if ~isempty(gw_data{2})
2697 soundData(2,:)=interp1(gw_data{2}(1,:),gw_data{2}(2,:)*opt.channelWeights(2),(0:(sound_length-1))*dwell_time,'linear',0);
2698 end
2699 if ~isempty(gw_data{3})
2700 tmp=interp1(gw_data{3}(1,:),0.5*gw_data{3}(2,:)*opt.channelWeights(3),(0:(sound_length-1))*dwell_time,'linear',0);
2701 soundData(1,:)=soundData(1,:)+tmp;
2702 soundData(2,:)=soundData(2,:)+tmp;
2703 end
2705 % filter like we did it in the gradient music project
2706 %b = fir1(40, 10000/sample_rate);
2707 %sound_data = filter(b, 1, sound_data,[],2);
2708 % use Gaussian convolution instead to supress ringing
2709 gw=gausswin(round(sample_rate/6000)*2+1);
2710 gw=gw/sum(gw(:));
2711 soundData(1,:) = conv(soundData(1,:), gw, 'same');
2712 soundData(2,:) = conv(soundData(2,:), gw, 'same');
2714 sound_data_max=max(abs(soundData(:)));
2715 soundData = 0.95 * soundData / sound_data_max;
2717 if ~opt.onlyProduceSoundData
2718 % info
2719 fprintf('playing out the sequence waveform, duration %.1gs\n', sound_length*dwell_time);
2720 % play out the sound
2721 % we have to zero-pad the weveform due to the limitations of
2722 % matlab-to-sound interface
2723 sound([zeros(2,sample_rate/2) soundData zeros(2,sample_rate/2)], sample_rate);
2724 end
2725 end
2727 function ok=install(seq,param1,param2)
2728 %install Install sequence on RANGE system.
2729 % install(seq) Install sequence by copying files to Siemens
2730 % host and RANGE controller
2731 %
2732 % install(seq,'sequence_path_or_name') Auto-detect scanner
2733 % environment and install the sequence under the given file
2734 % name. If sub-directories are provided prior to the name
2735 % they will be create automatically.
2736 %
2737 % install(seq,'siemens') Install Siemens Numaris4 file as external.seq
2738 % install(seq,'siemensNX') Install Siemens NumarisX file as external.seq
2739 % install(seq,'siemens','sequence_path_or_name') Install
2740 % Pulseq file assuming a Numaris4 Siemens system
2741 % under the given name and optinally path.
2742 % install(seq,'siemens','sequence_path_or_name') Install
2743 % Pulseq file assuming a NumarisX Siemens system
2744 % under the given name and optinally path.
2746 if ispc
2747 % windows
2748 ping_command='ping -w 1000 -n 1';
2749 elseif isunix || ismac
2750 % unix-like
2751 ping_command='ping -q -n -W1 -c1';
2752 end
2754 % for compatibility with older versions
2755 if nargin==3
2756 name=param2;
2757 dest=param1;
2758 else
2759 switch param1
2760 case {'siemens','siemensNX'...
2761 }
2762 name='external';
2763 dest=param1;
2764 otherwise
2765 name=param1;
2766 % auto-detect the scanner environment
2767 cmd=[ping_command ' 192.168.2.2 && ssh -oBatchMode=yes -oStrictHostKeyChecking=no -oHostKeyAlgorithms=+ssh-rsa root@192.168.2.2 ls /opt/medcom/MriCustomer/CustomerSeq'];
2768 [status, ~] = system(cmd);
2769 if status == 0
2770 fprintf('Siemens NumarisX environment detected\n');
2771 dest='siemensNX';
2772 else
2773 fprintf('Assuming Siemens Numaris4 environment (not tested yet)\n');
2774 dest='siemens';
2775 end
2776 end
2777 end
2779 ok = true;
2780 if any(strcmpi(dest, {'both', 'siemens', 'siemensNX'}))
2781 seq.write('external.seq.tmp');
2782 [filepath,filename,ext] = fileparts(name);
2783 filepath=strrep(filepath,'\','/'); % fix for windows users
2784 if contains(filepath,'../')
2785 error('No relative path elements (like ..) are allowed.');
2786 end
2787 name=[filepath '/' filename]; % discard the extension
2788 % we assume we are in the internal network but ICE computer
2789 % on VB and VE has different IPs
2790 if ~strcmpi(dest, 'siemensNX')
2791 ice_ips={'192.168.2.3', '192.168.2.2'};
2792 ice_ip=[];
2793 for i=1:length(ice_ips)
2794 [status, ~] = system([ping_command ' ' ice_ips{i}]);
2795 if status == 0
2796 ice_ip=ice_ips{i};
2797 break;
2798 end
2799 end
2800 else
2801 ice_ip='192.168.2.2';
2802 end
2803 if isempty(ice_ip)
2804 error('Scanner not found, sequence install failed.')
2805 end
2806 pulseq_seq_path='/opt/medcom/MriCustomer/seq/pulseq';
2807 if strcmpi(dest,'siemensNX')
2808 pulseq_seq_path='/opt/medcom/MriCustomer/CustomerSeq/pulseq';
2809 end
2810 [status, retmes] = system(['scp -oBatchMode=yes -oStrictHostKeyChecking=no -oHostKeyAlgorithms=+ssh-rsa external.seq.tmp root@' ice_ip ':' pulseq_seq_path '/external_tmp.seq']);
2811 ok = ok & status == 0;
2812 if ok
2813 if ~isempty(filepath)
2814 mkdir_add=['mkdir -p ' pulseq_seq_path '/' filepath ';'];
2815 else
2816 mkdir_add=[];
2817 end
2818 sys_cmd=['ssh -oBatchMode=yes -oStrictHostKeyChecking=no -oHostKeyAlgorithms=+ssh-rsa root@' ice_ip ' "chmod a+rw ' pulseq_seq_path '/external_tmp.seq ;' mkdir_add ' rm -f ' pulseq_seq_path '/' name '.seq; mv ' pulseq_seq_path '/external_tmp.seq ' pulseq_seq_path '/' name '.seq; ls -l ' pulseq_seq_path '/' name '.seq"'];
2819 fprintf('running command: %s\n',sys_cmd);
2820 [status,cmdout] = system(sys_cmd);
2821 %[status,cmdout] = system(['start cmd /c ssh -oBatchMode=yes -oStrictHostKeyChecking=no -oHostKeyAlgorithms=+ssh-rsa root@' ice_ip ' "rm -f ' pulseq_seq_path '/' name '.seq"']);
2822 %[status,cmdout] = system(['start cmd /c ssh -oBatchMode=yes -oStrictHostKeyChecking=no -oHostKeyAlgorithms=+ssh-rsa root@' ice_ip ' "mv ' pulseq_seq_path '/external_tmp.seq ' pulseq_seq_path '/' name '.seq"']);
2823 else
2824 error(['Failed to copy the sequence file to the scanner, the returned error message is: ' mr.aux.strstrip(retmes)]);
2825 end
2826 end
2829 if ok
2830 fprintf('Sequence installed as %s.seq\n',name)
2831 else
2832 error('Sequence install failed.')
2833 end
2834 end
2836 function id = getExtensionTypeID(obj, str)
2837 % get numeric ID for the given string extention ID
2838 % will automatically create a new ID if unknown
2839 num=find(strcmp(obj.extensionStringIDs,str));
2840 if isempty(num)
2841 if isempty(obj.extensionNumericIDs)
2842 id=1;
2843 else
2844 id=1+max(obj.extensionNumericIDs);
2845 end
2846 obj.extensionNumericIDs(1+length(obj.extensionNumericIDs))=id;
2847 obj.extensionStringIDs{1+length(obj.extensionStringIDs)}=str;
2848 assert(length(obj.extensionNumericIDs)==length(obj.extensionStringIDs));
2849 else
2850 id=obj.extensionNumericIDs(num);
2851 end
2852 end
2854 function str = getExtensionTypeString(obj, id)
2855 % get numeric ID for the given string extention ID
2856 % may fail
2857 num=find(obj.extensionNumericIDs==id);
2858 if isempty(num)
2859 error(['Extension for the given ID ' num2str(id) ' is unknown']);
2860 end
2861 str=obj.extensionStringIDs{num};
2862 end
2864 function setExtensionStringAndID(obj, str, id)
2865 % set numeric ID for the given string extention ID
2866 % may fail if not unique
2867 if any(strcmp(obj.extensionStringIDs,str)) || any(obj.extensionNumericIDs==id)
2868 error('Numeric or String ID is not unique');
2869 end
2870 obj.extensionNumericIDs(1+length(obj.extensionNumericIDs))=id;
2871 obj.extensionStringIDs{1+length(obj.extensionStringIDs)}=str;
2872 assert(length(obj.extensionNumericIDs)==length(obj.extensionStringIDs))
2873 end
2875 function id = getOrCreateTridId(obj, label_name)
2876 if isstring(label_name)
2877 label_name = char(label_name);
2878 end
2879 if ~ischar(label_name) || isempty(label_name)
2880 error('TRID label_name must be a non-empty char/string.');
2881 end
2883 if isKey(obj.tridName2Id, label_name)
2884 id = obj.tridName2Id(label_name);
2885 else
2886 id = int32(numel(obj.tridId2Name) + 1);
2887 obj.tridName2Id(label_name) = id;
2888 obj.tridId2Name{double(id),1} = label_name;
2889 end
2891 obj.tridHistory{end+1,1} = label_name;
2892 end
2893 end
2895 methods(Static)
2896 function codes = getBinaryCodes()
2897 %getBinaryCodes Return binary codes for section headers in
2898 % in a binary sequence file.
2899 %
2900 % See also writeBinary
2902 codes.fileHeader = typecast([uint8(1) uint8('pulseq') uint8(2)],'int64');
2903 %codes.version_major = int64(obj.version_major);
2904 %codes.version_minor = int64(obj.version_minor);
2905 %codes.version_revision = int64(obj.version_revision);
2906 prefix = bitshift(int64(hex2dec('FFFFFFFF')), 32);
2907 codes.section.definitions = bitor(prefix, int64(1));
2908 codes.section.blocks = bitor(prefix, int64(2));
2909 codes.section.rf = bitor(prefix, int64(3));
2910 codes.section.gradients = bitor(prefix, int64(4));
2911 codes.section.trapezoids = bitor(prefix, int64(5));
2912 codes.section.adc = bitor(prefix, int64(6));
2913 codes.section.delays = bitor(prefix, int64(7));
2914 codes.section.shapes = bitor(prefix, int64(8));
2915 codes.section.extensions = bitor(prefix, int64(9));
2916 codes.section.triggers = bitor(prefix, int64(10));
2917 codes.section.labelset = bitor(prefix, int64(11));
2918 codes.section.labelinc = bitor(prefix, int64(12));
2919 codes.section.softdelays = bitor(prefix, int64(13));
2920 codes.section.rfshims = bitor(prefix, int64(14));
2921 codes.section.rotations = bitor(prefix, int64(15));
2922 %
2923 codes.section.signature = bitor(prefix, int64(0x00FFFFFF));
2924 end
2925 end
2927end % classdef