concept-collection / seqlab
seqlab / src / engine / pulseq / +mr / restoreAdditionalShapeSamples.m
46 lines · 2.5 KBBlameHistoryRaw
1function [tt_chg, waveform_chg] = restoreAdditionalShapeSamples(tt,waveform,first,last,gradRasterTime,iBlock)
2% restore shape: if we had a
3% trapezoid converted to shape we have to find
4% the "corners" and we can eliminate internal
5% samples on the straight segments
6% but first we have to restore samples on the
7% edges of the gradient raster intervals
8% for that we need the first sample
9 max_abs=max(abs(waveform));
10 odd_step1=[first 2*waveform'];
11 odd_step2=odd_step1.*(mod(1:length(odd_step1),2)*2-1);
12 waveform_odd_rest=(cumsum(odd_step2).*(mod(1:length(odd_step2),2)*2-1))';
13 waveform_odd_interp=[first; 0.5*(waveform(1:end-1)+waveform(2:end)); last];
14 if abs(waveform_odd_rest(end)-last)>2e-5*max_abs % what's the reasonable threshold?
15 blInfo='';
16 if exist('iBlock')
17 blInfo=['[block ' num2str(iBlock) '] '];
18 end
19 warning('mr:restoreShape',[blInfo 'Last restored point ' ...
20 'differs too much from the recorded last, skipping the shape restoration step; ' ...
21 'deviation: ' num2str(abs(waveform_odd_rest(end)-last)) 'Hz/m (' num2str(abs(waveform_odd_rest(end)-last)/max_abs*100) '%%); ' ...
22 'No reason to panic, this is typical for spirals. ' ...
23 'To disable this warning (assuming you know what you are doing) add \n' ...
24 '`warning(''OFF'', ''mr:restoreShape'')ยด to your sequence-generating script.']);
25 tt_chg=[0 tt' tt(end)+gradRasterTime/2];
26 waveform_chg = [first waveform' last];
27 return;
28 end
29 %figure; plot([0,10e-6+grad.t'],waveform_odd_rest-waveform_odd_interp);
30 waveform_odd_mask=abs(waveform_odd_rest-waveform_odd_interp)<=eps+2e-5*max_abs; % threshold ???
31 waveform_odd=waveform_odd_interp.*waveform_odd_mask+waveform_odd_rest.*(1-waveform_odd_mask);
33 % combine odd & even
34 comb=[ 0 waveform' ; waveform_odd' ];
35 waveform_os=comb(2:end)';
37 tt_odd=(0:(length(waveform_odd_rest)-1))*gradRasterTime;
38 tt_os=(0:(length(waveform_os)-1))*gradRasterTime*0.5;
40 waveform_even_reint=0.5*(waveform_odd_rest(1:end-1)+waveform_odd_rest(2:end));
42 maskChanges = abs([1; diff(waveform_os,2); 1])>1e-8; % TRUE if values change
43 waveform_chg = waveform_os(maskChanges)'; % Elements without repetitions
44 tt_chg=tt_os(maskChanges);
45 %figure;plot(grad.tt,grad.waveform);hold on; plot(tt_chg,waveform_chg); plot(tt_chg,waveform_chg,'o');
46end