78d04f1seqlab: write and view pulseq MRI sequences in the browserJeremy Magland 1function [varargout] = splitGradientAt(grad, timepoint, varargin)
2%SplitGradient Splits a trapezoidal gradient into two extended trapezoids
3%(currently shaped gradients) defined by the cut line.
4%
5% [grads] = splitGradient(grad)
6% Returns the two gradient parts by cutting the original 'grad' at the
7% 'timepoint' . For the input type 'trapezoid' the results are trtyurned
8% as extended trapezoids, for 'arb' as arbitrary gradient objects. The
9% delays in the individual gradient events are adapted such that
10% addGradients(...) produces an gradient equivalent to 'grad'.
11%
12% See also splitGradient makeExtendedTrapezoid makeTrapezoid
13% Sequence.addBlock mr.opts
14%
15% Maxim Zaitsev <maxim.zaitsev@uniklinik-freiburg.de>
16% Stefan Kroboth <stefan.kroboth@uniklinik-freiburg.de>
18persistent parser
20if isempty(parser)
21 parser = mr.aux.InputParserCompat;
22 parser.FunctionName = 'splitGradientAt';
23 parser.addRequired('grad', @isstruct);
24 parser.addRequired('timepoint', @isnumeric);
25 parser.addOptional('system', [], @isstruct);
26end
27parse(parser, grad, timepoint, varargin{:});
28opt = parser.Results;
30if isfield(grad,'id')
31 error('attempting to split readily registered object! please register objects after calling this function or deregister the argument by calling rmfield(...,''id'')');
32end
34if isempty(opt.system)
35 system=mr.opts();
36else
37 system=opt.system;
38end
40gradRasterTime = system.gradRasterTime;
42% round the time point to the gradient raster;
43timeindex = round(timepoint / gradRasterTime);
44if abs(timepoint-timeindex*gradRasterTime)>1e-6
45 warning('splitting the gradint at a point that is not on a gradient raster edge, substantial rounding is applied');
46end
47timepoint = timeindex * gradRasterTime;
48timeindex = timeindex + 1; % convert to Matlab convention
50ch = grad.channel;
52if strcmp(grad.type, 'grad')
53 % check if we have an arbitrary gradient or an exended trapezoid
54 if abs(grad.tt(1)-0.5*gradRasterTime)<1e-10
55 % it can be an arbitrary gradient or arbitrary gradient with oversampling
56 isArb=all(abs(grad.tt(2:end)-grad.tt(1:end-1)-gradRasterTime)<1e-10);
57 isArbOs=all(abs(grad.tt(2:end)-grad.tt(1:end-1)-gradRasterTime*0.5)<1e-10);
58 if isArb || isArbOs
59 if isArbOs
60 % update timeindex to account for dencier sampling
61 timeindex = (timeindex-1)*2;
62 end
63 % arbitrary gradient -- the most trivial conversion
64 % if timepoint is out of range we have nothing to do
65 if timeindex == 1 || timeindex >= length(grad.tt)
66 varargout{1} = grad;
67 else
68 grad1=grad;
69 grad2=grad;
70 if isArbOs
71 grad1.last=grad.waveform(timeindex);
72 else
73 grad1.last=0.5*(grad.waveform(timeindex-1)+grad.waveform(timeindex)); % FIXME: retrive the double-sampling point (e.g. the corner of the trapezoid)
74 end
75 grad2.first=grad1.last;
76 grad2.delay=grad.delay + timepoint;
77 grad1.tt=grad.tt(1:(timeindex-1));
78 grad1.waveform=grad.waveform(1:(timeindex-1));
79 if isArbOs
80 grad2.tt=grad.tt(timeindex+1:end) - timepoint;
81 grad2.waveform=grad.waveform(timeindex+1:end);
82 else
83 grad2.tt=grad.tt(timeindex:end) - timepoint;
84 grad2.waveform=grad.waveform(timeindex:end);
85 end
86 grad1.shape_dur = grad1.tt(end) - grad1.tt(1) + gradRasterTime;
87 grad2.shape_dur = grad2.tt(end) - grad2.tt(1) + gradRasterTime;
89 if nargout==1
90 varargout{1} = [grad1 grad2];
91 else
92 varargout{1} = grad1;
93 varargout{2} = grad2;
94 end
95 end
96 %figure; plot(grad.tt, grad.waveform); hold on; plot(grad1.tt, grad1.waveform); plot(grad2.delay+grad2.tt, grad2.waveform);
97 return; % early return to protect the subsequent code
98 end
99 end
101 % we have an extended trapezoid (by excluding arbitrary grad) -- excellent choice!
102 times = grad.tt';
103 amplitudes = grad.waveform'; % QC: to match the matrix size for times1 and amplitudes1 below. 2025.01.02
105elseif strcmp(grad.type, 'trap')
106 grad.delay = round(grad.delay /gradRasterTime)*gradRasterTime; % MZ: was ceil
107 grad.riseTime = round(grad.riseTime/gradRasterTime)*gradRasterTime; % MZ: was ceil
108 grad.flatTime = round(grad.flatTime/gradRasterTime)*gradRasterTime; % MZ: was ceil
109 grad.fallTime = round(grad.fallTime/gradRasterTime)*gradRasterTime; % MZ: was ceil
111 % prepare the extended trapezoid structure
112 if grad.flatTime == 0
113 times = [0 grad.riseTime grad.riseTime+grad.fallTime];
114 amplitudes = [0 grad.amplitude 0];
115 else
116 times = [0 grad.riseTime grad.riseTime+grad.flatTime grad.riseTime+grad.flatTime+grad.fallTime];
117 amplitudes = [0 grad.amplitude grad.amplitude 0];
118 end
119else
120 error('Splitting of unsupported event.');
121end
123% if the cutline is behind the gradient there is no second gradient to create
124if timepoint >= grad.delay+times(end)
125 error('trying to place the splitting time point after the end of the gradient');
126end
128% now we have everything in the extended trapezoid structure
130% if the cutline goes through the delay we need special treatment
131if timepoint < grad.delay
132 times=[0 grad.delay+times];
133 amplitudes = [0 amplitudes];
134 grad.delay=0;
135else
136 timepoint = timepoint - grad.delay;
137end
139% sample at timepoint
140amp_tp=interp1(times, amplitudes, timepoint, 'linear'); % MZ: interp1() is not OK here for the corner situation TODO: fixme! (e.g. by restoring the corners as done in waveforms_and_times())
141% split the data
142teps=1e-10; % we need this because of the rounding problems
143times1 = [ times(times<timepoint-teps) timepoint ];
144amplitudes1 = [ amplitudes(times<timepoint-teps) amp_tp ];
145times2 = [ timepoint times(times>timepoint+teps) ] - timepoint;
146amplitudes2 = [ amp_tp amplitudes(times>timepoint+teps) ];
148% recreate gradients
149grad1 = mr.makeExtendedTrapezoid(ch, 'system', system, 'times', times1,...
150 'amplitudes', amplitudes1, ...
151 'skip_check', true);
152grad1.delay = grad.delay;
153grad2 = mr.makeExtendedTrapezoid(ch, 'system', system, 'times', times2,...
154 'amplitudes', amplitudes2, ...
155 'skip_check', true);
156grad2.delay = timepoint + grad.delay;
158%grads = [grad1 grad2];
159if nargout==1
160 varargout{1} = [grad1 grad2];
161else
162 varargout{1} = grad1;
163 varargout{2} = grad2;
164end
166end