78d04f1seqlab: write and view pulseq MRI sequences in the browserJeremy Magland 1classdef SeqPlot < handle
2 %plot Plot the sequence in a new figure.
3 % plot(seqObj) Plot the sequence
4 %
5 % plot(...,'timeRange',[start stop]) Plot the sequence
6 % between the times specified by start and stop.
7 %
8 % plot(...,'blockRange',[first last]) Plot the sequence
9 % starting from the first specified block to the last one.
10 %
11 % plot(...,'timeDisp',unit) Display time in:
12 % 's', 'ms' or 'us'.
13 %
14 % plot(...,'label','LIN,REP') Plot label values for ADC events:
15 % in this example for LIN and REP labels; other valid labes are
16 % accepted as a comma-separated list.
17 %
18 % plot(...,'showBlocks',1) Plot grid and tick labels at the
19 % block boundaries. Accepts a numeric or a boolean parameter.
20 %
21 % plot(...,'stacked',1) Rearrange the plots such they are vertically
22 % stacked and share the same x-axis. Accepts a numeric or a boolean
23 % parameter.
24 %
25 % plot(...,'showGuides',1) How dynamic hairline guides that follow
26 % the data cursor to help verifying event alignment. Accepts a
27 % numeric or a boolean parameter.
28 %
29 % f=plot(...) Return the new figure handle.
30 %
32 properties (Access = public)
33 f % figure handle
34 end
36 properties (Access = private)
37 ax % array of plot axes handles
38 vLines % array of vline handles
40 hSeq
42 end
44 properties (Constant = true, Hidden = true)
46 % vertical margin (px)
47 margin = 6;
48 % lower vertical margin (px)
49 my1 = 45;
50 % left horizontal margin
51 mx1 = 70;
52 % right horizontal margin
53 mx2 = 5;
55 end
58 methods
60 function obj = SeqPlot(seq, varargin)
62 validTimeUnits = {'s','ms','us'};
63 validLabel = mr.getSupportedLabels();
64 persistent parser
65 if isempty(parser)
66 parser = inputParser;
67 parser.FunctionName = 'plot';
68 parser.addParamValue('showBlocks',false,@(x)(isnumeric(x) || islogical(x)));
69 parser.addParamValue('timeRange',[0 inf],@(x)(isnumeric(x) && length(x)==2));
70 parser.addParamValue('blockRange',[1 inf],@(x)(isnumeric(x) && length(x)==2));
71 parser.addParamValue('timeDisp',validTimeUnits{1},...
72 @(x) any(validatestring(x,validTimeUnits)));
73 parser.addParamValue('label',[]);%,@(x)(isstr(x)));%@(x) any(validatestring(x,validLabel))
74 parser.addParamValue('hide',false);%,@(x)(isstr(x)));%@(x) any(validatestring(x,validLabel))
75 parser.addParamValue('stacked',false);%,@(x)(isstr(x)));%@(x) any(validatestring(x,validLabel))
76 parser.addParamValue('showGuides',true);%,@(x)(isstr(x)));%@(x) any(validatestring(x,validLabel))
77 end
78 parse(parser,varargin{:});
79 opt = parser.Results;
81 if mr.aux.isOctave()
82 if opt.stacked
83 warning('Option stacked is not (yet) supported by Octave');
84 opt.stacked=false;
85 end
86 if opt.showBlocks
87 warning('Option stacked is not (yet) supported by Octave');
88 opt.showBlocks=false;
89 end
90 end
92 obj.f=figure;
93 obj.hSeq=seq; % Sequence is a handle-class so copying is cheap...
95 set(obj.f, 'Visible', 'off')
97 if ~mr.aux.isOctave()
98 obj.ax = gobjects(1,6);
99 end
100 for i=1:6
101 obj.ax(i)=subplot(3,2,i);
102 end
103 obj.ax=obj.ax([1 3 5 2 4 6]); % Re-order axes
104 arrayfun(@(x)hold(x,'on'),obj.ax);
105 arrayfun(@(x)grid(x,'on'),obj.ax);
106 labels={'ADC/lbl/trig','RF mag (Hz)','RF/ADC ph (rad)','Gx (kHz/m)','Gy (kHz/m)','Gz (kHz/m)'};
107 arrayfun(@(x)ylabel(obj.ax(x),labels{x}),1:6);
109 tFactorList = [1 1e3 1e6];
110 tFactor = tFactorList(strcmp(opt.timeDisp,validTimeUnits));
111 xlabel(obj.ax(3),['t (' opt.timeDisp ')']);
112 xlabel(obj.ax(6),['t (' opt.timeDisp ')']);
114 t0=0;
115 label_defined=false;
116 label_indexes_2plot=[];
117 label_legend_2plot=[];
118 for i=1:length(validLabel)
119 label_store.(validLabel{i})=0;
120 if ~isempty(opt.label) && ~isempty(strfind(upper(opt.label),validLabel{i}))
121 label_indexes_2plot=[label_indexes_2plot i];
122 label_legend_2plot=[label_legend_2plot; validLabel{i}];
123 end
124 end
125 if ~isempty(label_indexes_2plot)
126 if mr.aux.isOctave()
127 label_colors_2plot=turbo(length(label_indexes_2plot)+1); % need +1 because the ADC plot by itself also "eats up" one color
128 else
129 label_colors_2plot=parula(length(label_indexes_2plot)+1); % need +1 because the ADC plot by itself also "eats up" one color
130 end
131 label_colors_2plot=[label_colors_2plot(end,:); label_colors_2plot(1:end-1,:)]; % we like these colors better ?
132 end
134 % time format
135 switch opt.timeDisp
136 case 'us'
137 timeFormat='%.1f';
138 case 'ms'
139 timeFormat='%.4f';
140 otherwise
141 timeFormat='%.7f';
142 end
144 % data cursor callback
145 if ~mr.aux.isOctave()
146 hDCM = datacursormode(obj.f);
147 hDCM.UpdateFcn = @(src, event)DataTipHandler(obj,tFactor,[timeFormat ' ' opt.timeDisp],src,event);
148 end
150 % time/block range
151 timeRange=opt.timeRange;
152 blockEdges=[0 cumsum(seq.blockDurations)];
153 if opt.blockRange(1)>1 && blockEdges(opt.blockRange(1))>timeRange(1)
154 timeRange(1)=blockEdges(opt.blockRange(1));
155 end
156 if isfinite(opt.blockRange(2)) && opt.blockRange(2)<length(seq.blockDurations) && blockEdges(opt.blockRange(2)+1)<timeRange(2)
157 timeRange(2)=blockEdges(opt.blockRange(2)+1);
158 end
159 % block timings
160 blockEdgesInRange=blockEdges(logical((blockEdges>=timeRange(1)).*(blockEdges<=timeRange(2))));
161 if strcmp(opt.timeDisp,'us') && ~mr.aux.isOctave()
162 for i=1:6
163 xax=get(obj.ax(i),'XAxis');
164 xax.ExponentMode='manual';
165 xax.Exponent=0;
166 end
167 end
168 if opt.showBlocks
169 % show block edges in plots
170 for i=1:6
171 xax=get(obj.ax(i),'XAxis');
172 xax.TickValues=unique(tFactor.*blockEdgesInRange);
173 set(obj.ax(i),'XTickLabelRotation',90);
174 %xax.MinorTickValues=tFactor.*blockEdgesInRange;
175 %set(obj.ax(i),'XMinorTick', 'on');
176 %set(obj.ax(i),'XMinorGrid', 'on');
177 %set(obj.ax(i),'GridColor',0.8*[1 1 1]);
178 %set(obj.ax(i),'MinorGridColor',0.6*[1 1 1]);
179 %set(obj.ax(i),'MinorGridLineStyle','-');
180 end
181 end
182 %
183 gradChannels={'gx','gy','gz'};
185 % loop through blocks
186 for iB=1:length(seq.blockEvents)
187 block = seq.getBlock(iB);
188 if isfield(block,'rotation')
189 % apply the rotation to the current block and restore the block structure
190 c=mr.rotate3D(block.rotation.rotQuaternion,block,'system',seq.sys);
191 for i=1:3
192 block.(gradChannels{i})=[];
193 end
194 for i=1:length(c)
195 if isstruct(c{i}) && isfield(c{i},'type') && isfield(c{i},'channel')
196 block.(['g' c{i}.channel])=c{i};
197 end
198 end
199 end
200 if t0<=timeRange(2)
201 % update the labels / counters even if we are below the display range
202 if isfield(block,'label') %current labels, works on the curent or next adc
203 for i=1:length(block.label)
204 if strcmp(block.label(i).type,'labelinc')
205 label_store.(block.label(i).label)=...
206 label_store.(block.label(i).label)+block.label(i).value;
207 else
208 label_store.(block.label(i).label)=block.label(i).value;
209 end
210 end
211 label_defined=true;
212 end
213 end
214 isValid = t0+seq.blockDurations(iB)>timeRange(1) && t0<=timeRange(2);
215 if isValid
216 if isfield(block,'trig') && ~isempty(block.trig)
217 switch(block.trig.type)
218 case 'output'
219 % plot digital output triggers in the RF-TX pane
220 p2x=plot(tFactor*(t0+block.trig.delay),0,'diamond','Color',[0 0.5 0],'Parent',obj.ax(1));
221 p2x=plot(tFactor*(t0+block.trig.delay +[0 block.trig.duration]),[0 0],'-','Marker','.','Color',[0 0.5 0],'Parent',obj.ax(1));
222 case 'trigger'
223 p1x=plot(tFactor*(t0+block.trig.delay),0,'>b','Parent',obj.ax(1));
224 p1x=plot(tFactor*(t0+block.trig.delay),0,'.b','Parent',obj.ax(1));
225 %otherwise
226 end
227 end
228 if ~isempty(block.adc)
229 adc=block.adc;
230 t=adc.delay + ((0:adc.numSamples-1)'+0.5)*adc.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)
231 p1=plot(tFactor*(t0+t),zeros(size(t)),'rx','Parent',obj.ax(1));
232 if isempty(adc.phaseModulation)
233 adc.phaseModulation=0;
234 end
235 full_freqOffset=adc.freqOffset+adc.freqPPM*1e-6*seq.sys.gamma*seq.sys.B0;
236 full_phaseOffset=adc.phaseOffset+adc.phasePPM*1e-6*seq.sys.gamma*seq.sys.B0;
237 p2=plot(tFactor*(t0+t), angle(exp(1i*(full_phaseOffset+adc.phaseModulation)).*exp(1i*2*pi*t*full_freqOffset)),'b.','MarkerSize',1,'Parent',obj.ax(3)); % plot ADC phase
238 % labels/counters/flags
239 if label_defined && ~isempty(label_indexes_2plot)
240 set(obj.ax(1),'ColorOrder',label_colors_2plot);
241 label_store_cell=struct2cell(label_store);
242 lbl_vals=[label_store_cell{label_indexes_2plot}];
243 t=t0+adc.delay + (adc.numSamples-1)/2*adc.dwell;
244 p=plot(tFactor*t,lbl_vals,'.','markersize',5,'Parent',obj.ax(1));
245 if ~isempty(label_legend_2plot)
246 legend(obj.ax(1),p,label_legend_2plot,'location','Northwest','AutoUpdate','off');
247 label_legend_2plot=[];
248 end
249 end
250 end
251 if ~isempty(block.rf)
252 rf=block.rf;
253 [tc,ic,fi]=mr.calcRfCenter(rf);
254 if fi==0
255 sc=rf.signal(ic);
256 else
257 sc=rf.signal(ic)*(1-abs(fi))+rf.signal(ic+sign(fi))*abs(fi);
258 end
259 if max(abs(diff(rf.t)-rf.t(2)+rf.t(1)))<1e-9 && length(rf.t)>100
260 % homogeneous sampling and long pulses -- use lower time resolution for better display and performance
261 dt=rf.t(2)-rf.t(1);
262 st=max(1,round(seq.sys.gradRasterTime/dt));
263 t=rf.t(1:st:end);
264 s=rf.signal(1:st:end);
265 % always include the last point for the accurate display
266 if (t(end)~=rf.t(end))
267 t(end+1)=rf.t(end);
268 s(end+1)=rf.signal(end);
269 end
270 else
271 t=rf.t;
272 s=rf.signal;
273 end
274 sreal=max(abs(imag(s)))/max(abs(real(s)))<1e-6; %all(isreal(s));
275 full_freqOffset=rf.freqOffset+rf.freqPPM*1e-6*seq.sys.gamma*seq.sys.B0;
276 full_phaseOffset=rf.phaseOffset+rf.phasePPM*1e-6*seq.sys.gamma*seq.sys.B0;
277 % If off-resonant and rectangular (2 samples), interpolate the pulse
278 if (length(s) == 2) && (full_freqOffset ~= 0)
279 numInterp = min(int32(abs(full_freqOffset)), 256);
280 t = linspace(t(1), t(end), numInterp)';
281 s = linspace(s(1), s(end), numInterp)';
282 end
283 if abs(s(1))~=0 % fix strangely looking phase / amplitude in the beginning
284 s=[0; s];
285 t=[t(1); t];
286 %ic=ic+1;
287 end
288 if abs(s(end))~=0 % fix strangely looking phase / amplitude at the end
289 s=[s; 0];
290 t=[t; t(end)];
291 end
293 if (sreal)
294 p1=plot(tFactor*(t0+t+rf.delay), real(s),'Parent',obj.ax(2));
295 p2=plot(tFactor*(t0+t+rf.delay), angle(s.*sign(real(s))*exp(1i*full_phaseOffset).*exp(1i*2*pi*t *full_freqOffset)), tFactor*(t0+tc+rf.delay), angle(sc*exp(1i*full_phaseOffset).*exp(1i*2*pi*tc*full_freqOffset)),'xb', 'Parent',obj.ax(3));
296 else
297 p1=plot(tFactor*(t0+t+rf.delay), abs(s),'Parent',obj.ax(2));
298 p2=plot(tFactor*(t0+t+rf.delay), angle(s*exp(1i*full_phaseOffset).*exp(1i*2*pi*t *full_freqOffset)), tFactor*(t0+tc+rf.delay), angle(sc*exp(1i*full_phaseOffset).*exp(1i*2*pi*tc*full_freqOffset)),'xb', 'Parent',obj.ax(3));
299 end
300 end
301 for j=1:length(gradChannels)
302 grad=block.(gradChannels{j});
303 if ~isempty(grad)
304 if strcmp(grad.type,'grad')
305 % we extend the shape by adding the first
306 % and the last points in an effort of
307 % making the display a bit less confusing...
308 %t=grad.delay + [0; grad.t + (grad.t(2)-grad.t(1))/2; grad.t(end) + grad.t(2)-grad.t(1)];
309 t= grad.delay+[0; grad.tt; grad.shape_dur];
310 waveform=1e-3* [grad.first; grad.waveform; grad.last];
311 else
312 t=cumsum([0 grad.delay grad.riseTime grad.flatTime grad.fallTime]);
313 waveform=1e-3*grad.amplitude*[0 0 1 1 0];
314 end
315 p=plot(tFactor*(t0+t),waveform,'Parent',obj.ax(3+j));
316 end
317 end
318 end
319 t0=t0+seq.blockDurations(iB);%mr.calcDuration(block);
320 end
322 % Set axis limits and zoom properties
323 dispRange = tFactor*[timeRange(1) min(timeRange(2),t0)];
324 arrayfun(@(x)xlim(x,dispRange),obj.ax);
325 linkaxes(obj.ax(:),'x')
326 if ~mr.aux.isOctave()
327 h = zoom(obj.f);
328 setAxesZoomMotion(h,obj.ax(1),'horizontal');
329 end
330 % manually fix the phase vertical scale to +- pi
331 ylim(obj.ax(3),[-pi pi]);
332 % make Y-axes little bit less tight
333 arrayfun(@(x) ylim(x, ylim(x) + 0.03*[-1 1]*sum(ylim(x).*[-1 1])), obj.ax(2:end));
336 if opt.showGuides
337 if mr.aux.isOctave()
338 warning('Option showGuides is not implemented in Octave');
339 else
340 % add vertical lines and make them follow the cursor
341 % x-position
342 for ii = 1:numel(obj.ax)
343 obj.vLines(ii) = xline(obj.ax(ii), 0, 'r--');
344 end
345 end
346 end
348 if opt.stacked
349 % vertical stacking is defined in guiResize
350 set(obj.f, 'ResizeFcn', @obj.guiResize)
351 obj.guiResize()
352 end
354 if ~opt.hide
355 set(obj.f, 'Visible', 'on')
356 end
358 % do not assign to 'ans' when called without assigned variable
359 if nargout == 0
360 clear obj
361 end
362 end
364 function guiResize(obj, ~, ~)
365 % guiResize()
366 % Is called whenever the figure-shape is changed and makes
367 % sure all UI elements are correctly psotitioned. This
368 % function implements a vertical stacking of the individual
369 % axes.
371 nAxes = numel(obj.ax);
372 width = obj.f.Position(3);
373 height = obj.f.Position(4);
375 axHeight = (height - (nAxes-1)*obj.margin - obj.my1) / nAxes;
376 axWidth = width - obj.mx1 - obj.mx2;
378 for ii = 1:nAxes
379 set(obj.ax(ii), 'units', 'pixels', 'Position', [obj.mx1, height-ii*axHeight-(ii-1)*obj.margin, axWidth, axHeight])
380 if ii ~= nAxes
381 set(obj.ax(ii), 'Xlabel', [])
382 set(obj.ax(ii), 'XTickLabel', {})
383 end
384 end
385 end
387 function out=DataTipHandler(obj, tfactor, timeFormat, src, event)
388 if ~isa(event,'matlab.graphics.internal.DataTipEvent') || ...
389 ~isprop(event, 'Position') || length(event.Position)<2 || ...
390 ~isprop(event, 'Target')
391 out=[];
392 return;
393 end
394 ax=src.Host.Parent;
395 % get the relevant target from the y-axes title
396 at=lower(ax.YLabel.String);
397 if strcmp(at(1:3),'adc') || ...
398 (strcmp(at(1:6),'rf/adc') && strcmp(event.Target.LineStyle,'none') && strcmp(event.Target.Marker,'.')) % we need to check whether we are dealing with the ADC phase, which is also shown in the same panel as the RF
399 field='adc';
400 else
401 field=at(1:2);
402 end
403 % create the custom data tip as tex-formatted cell array of lines
404 t=event.Position(1);
405 t0=t;
406 if isa(event.Target,'matlab.graphics.chart.primitive.Line')
407 % for trapezoid gradients the last point may belong to the next block
408 t0=event.Target.XData(1);
409 end
410 iB=obj.hSeq.findBlockByTime(t0/tfactor);
411 rb=obj.hSeq.getRawBlockContentIDs(iB);
412 out={['\bf\color{blue}t:\rm\color{black}' sprintf(timeFormat,t)],...
413 ['\bf\color{blue}Y:\rm\color{black}' num2str(event.Position(2))],...
414 ''};
415 if isempty(rb.(field))
416 out{3}=['\bf\color{blue}blk:\rm\color{black}' num2str(iB)];
417 % we could add handling of the trigger/label data tips here
418 % specifically for the adc panel but it would imply a
419 % substantial performance hit because we'd have to unpack
420 % extensions, etc...
421 else
422 try
423 switch field(1)
424 case 'a'
425 name = obj.hSeq.adcID2NameMap(rb.(field));
426 case 'r'
427 name = obj.hSeq.rfID2NameMap(rb.(field));
428 otherwise
429 name = obj.hSeq.gradID2NameMap(rb.(field));
430 end
431 out{3}=['\bf\color{blue}blk:\rm\color{black}' num2str(iB) ' \bf\color{blue}' field '\_id:\rm\color{black}' num2str(rb.(field)) ' ''\bf\color{darkGreen}' name '\rm\color{black}'''];
432 catch
433 out{3}=['\bf\color{blue}blk:\rm\color{black}' num2str(iB) ' \bf\color{blue}' field '\_id:\rm\color{black}' num2str(rb.(field))];
434 end
435 end
437 % we need to delay the call of the update, otherwise the plot
438 % object generates an exception
439 t=timer('StartDelay',0e-3,'Period',1e-3,'TimerFcn',@(~,~)updateGuides(obj,t));
440 t.start();
441 end
443 function updateGuides(obj, tPos)
444 % updateGuides(tPos)
445 % updates the time-position for all vertical line objects in
446 % all axes
448 for ii = 1:numel(obj.vLines)
449 set(obj.vLines(ii), 'Value', tPos);
450 end
451 end
452 end
453end