1function grad = addGradients(grads, varargin)
2%addGradients Superposition of several gradients on a single axis.
3%
4% PURPOSE
5% Combines two or more gradient events on the same logical channel
6% ('x', 'y', or 'z') into a single equivalent gradient event by
7% pointwise summation of their waveforms on a common time grid. Used
8% to merge pre-phasers, spoilers, blip-up/blip-down pairs, and similar
9% co-scheduled gradients into one block entry consumed by
10% mr.Sequence/addBlock.
11%
12% The returned struct's type depends on the inputs:
13% - if every input is a trapezoid with identical delay, riseTime,
14% flatTime, and fallTime, the result is a trapezoid struct with
15% amplitudes summed (fast path);
16% - if every input is a trapezoid or an extended trapezoid on an
17% irregular time sampling, the result is an extended trapezoid
18% built via mr.makeExtendedTrapezoid;
19% - otherwise the result is an arbitrary gradient on the gradient
20% raster built via mr.makeArbitraryGrad.
21% In all three cases the returned struct is a valid gradient event
22% and can be passed directly to mr.Sequence/addBlock.
23%
24% SIGNATURES
25% grad = mr.addGradients(grads) % uses defaults from mr.opts()
26% grad = mr.addGradients(grads, system) % positional system
27% grad = mr.addGradients(grads, 'system', system) % name-value system
28% grad = mr.addGradients(grads, system, 'maxGrad', g, ...) % override slew/amplitude caps
29%
30% grads must be a cell array of at least two gradient event structs,
31% all on the same channel. Per-gradient delays are preserved; the
32% returned gradient's delay is the smallest delay among the inputs.
33%
34% INPUTS
35% grads [required] cell array of >=2 gradient event structs on the same channel
36% system [optional] struct from mr.opts. If omitted or empty, mr.opts() defaults
37% are used. Also accepted as 'system', sys name/value pair.
38% maxGrad [name/value] double, Hz/m, override for max gradient amplitude used by
39% the arbitrary-grad path. Default 0 (use system.maxGrad).
40% maxSlew [name/value] double, Hz/m/s, override for max slew rate used by the
41% arbitrary-grad path. Default 0 (use system.maxSlew).
42%
43% OUTPUT
44% grad struct. Gradient yype depends on the input mix (see PURPOSE).
45% Field order (this detailed description will be moved to respective mr.make... functions):
46%
47% Trapezoid fast path (all inputs identical-timing traps):
48% .type 'trap'
49% .channel 'x' | 'y' | 'z'
50% .amplitude Hz/m, sum of input amplitudes
51% .riseTime seconds, same as inputs
52% .flatTime seconds, same as inputs
53% .fallTime seconds, same as inputs
54% .area 1/m, sum of input areas
55% .flatArea 1/m, sum of input flat areas
56% .delay seconds, same as inputs
57% .first 0 (trapezoids always start at zero)
58% .last 0 (trapezoids always end at zero)
59%
60% Extended trapezoid path (all trap or extended-trap inputs):
61% .type 'grad'
62% .channel 'x' | 'y' | 'z'
63% .waveform Hz/m, amplitude samples on the union time grid
64% .delay seconds, min delay across inputs
65% .tt seconds, time-offsets of samples relative to delay
66% .shape_dur seconds, waveform duration on the gradient raster
67% .area 1/m, waveform area
68% .first Hz/m, sum of inputs' first values that share the common delay
69% .last Hz/m, sum of inputs' last values that share the max duration
70%
71% Arbitrary gradient fallback (mixed arbitrary + other):
72% .type 'grad'
73% .channel 'x' | 'y' | 'z'
74% .waveform Hz/m, uniformly-sampled amplitudes on system.gradRasterTime
75% (or half that if any input is oversampled-arbitrary)
76% .delay seconds, min delay across inputs
77% .area 1/m, waveform area
78% .tt seconds, sample time-offsets relative to delay
79% .shape_dur seconds, waveform duration
80% .first Hz/m, as above
81% .last Hz/m, as above
82%
83% ERRORS
84% - 'gradients have to be passed as cell array': grads is not a cell
85% - 'cannot add less then two gradients': numel(grads) < 2
86% - 'cannot add gradients on different channels': inputs mix x/y/z
87% Additional errors may propagate from mr.makeArbitraryGrad (e.g., slew
88% rate or amplitude limit violation) when the arbitrary-grad fallback
89% is taken.
90%
91% NOTES
92% - 'system' is registered as a positional (addOptional) parameter but
93% Pulseq's permissive input parser also accepts it as name-value.
94% - The returned delay is the minimum delay among the inputs; shapes
95% that start later than that are zero-padded internally so the
96% summed waveform preserves each input's original onset time.
97% - The fast trapezoid path triggers only when all inputs share the
98% same delay, riseTime, flatTime, AND fallTime. Inputs with matching
99% total duration but different rise/fall times (typical when the
100% caller uses 'Duration' + different 'Area') fall through to the
101% extended-trapezoid path.
102% - If any input is an oversampled arbitrary gradient, the result is
103% sampled at system.gradRasterTime/2 instead of system.gradRasterTime.
104% - maxGrad / maxSlew are only consulted on the arbitrary-grad path;
105% they are ignored on the trap and extended-trap paths (may change in future).
106%
107% EXAMPLE
108% sys = mr.opts('MaxGrad', 30, 'GradUnit', 'mT/m', ...
109% 'MaxSlew', 170, 'SlewUnit', 'T/m/s');
110%
111% % pre-phaser immediately followed 5 ms later by a spoiler, combined
112% % into one readout-axis gradient event
113% gxPre = mr.makeTrapezoid('x', 'Area', -500, 'Duration', 2e-3, 'system', sys);
114% gxSpoil = mr.makeTrapezoid('x', 'Area', 2000, 'Duration', 2e-3, ...
115% 'delay', 5e-3, 'system', sys);
116% gxComb = mr.addGradients({gxPre, gxSpoil}, 'system', sys);
117%
118% seq = mr.Sequence(sys);
119% seq.addBlock(gxComb);
120%
121% SEE ALSO
122% mr.Sequence/addBlock, mr.opts, mr.makeTrapezoid,
123% mr.makeExtendedTrapezoid, mr.makeArbitraryGrad, mr.calcDuration
124%
125% Maxim Zaitsev <maxim.zaitsev@uniklinik-freiburg.de>
126% Stefan Kroboth <stefan.kroboth@uniklinik-freiburg.de>
128persistent parser
130if isempty(parser)
131 parser = mr.aux.InputParserCompat;
132 parser.FunctionName = 'addGradients';
133 parser.addRequired('grads');
134 parser.addOptional('system', [], @isstruct);
135 parser.addParamValue('maxGrad', 0, @isnumeric);
136 parser.addParamValue('maxSlew', 0, @isnumeric);
137end
138parse(parser, grads, varargin{:});
139opt = parser.Results;
141if isempty(opt.system)
142 system=mr.opts();
143else
144 system=opt.system;
145end
147maxSlew = system.maxSlew;
148maxGrad = system.maxGrad;
149if opt.maxGrad > 0
150 maxGrad = opt.maxGrad;
151end
152if opt.maxSlew > 0
153 maxSlew = opt.maxSlew;
154end
156if ~iscell(grads)
157 error('gradients have to be passed as cell array');
158end
160if length(grads)<2
161 error('cannot add less then two gradients');
162end
164% first gradient event defines channel:
165channel = grads{1}.channel;
167% find out the general delay of all gradients and other statistics
168delays = []; % TODO: preallocate instead of grow
169firsts = [];
170lasts = [];
171durs=[];
172is_trap=[];
173is_arb=[];
174is_osa=[]; % oversampled arbitrary grad
175for ii = 1:length(grads)
176 if grads{ii}.channel~=channel
177 error('cannot add gradients on different channels');
178 end
179 delays = [delays, grads{ii}.delay];
180 durs = [durs, mr.calcDuration(grads{ii})];
181 is_trap = [is_trap, strcmp(grads{ii}.type,'trap')];
182 if is_trap(end)
183 is_arb = [is_arb, false];
184 is_osa = [is_osa, false];
185 % remember first/last
186 firsts = [firsts, 0];
187 lasts = [lasts, 0];
188 else
189 % check if this is an extended trapezoid
190 tt_rast=grads{ii}.tt/system.gradRasterTime;
191 is_arb = [is_arb, all(abs(tt_rast(:)'+0.5-(1:length(tt_rast)))<1e-6)];
192 is_osa = [is_osa, all(abs(tt_rast(:)'-0.5*(1:length(tt_rast)))<1e-6)];
193 % remember first/last
194 firsts = [firsts, grads{ii}.first];
195 lasts = [lasts, grads{ii}.last];
196 end
197end
198common_delay = min(delays);
199total_duration = max(durs);
200is_etrap=(~is_trap)&(~is_arb)&(~is_osa);
202% check if we have a set of traps with the same timing
203if all(is_trap)
204 % now all fields are the same so we can convert cell to a normal array
205 gradsa=cell2mat(grads);
206 if 1==length(unique([gradsa.delay])) && ...
207 1==length(unique([gradsa.riseTime])) && ...
208 1==length(unique([gradsa.flatTime])) && ...
209 1==length(unique([gradsa.fallTime]))
210 % TADA, all our gradients have the same timing, so we just add
211 % the amplitudes!
212 grad=gradsa(1);
213 grad.amplitude = sum([gradsa.amplitude]);
214 grad.area = sum([gradsa.area]);
215 grad.flatArea = sum([gradsa.flatArea]);
216 return;
217 end
218end
220% check if we only have arbitrary grads on irregular time samplings
221% optionally mixed with trapezoids
222if all(is_trap | is_etrap)
223 % we can do quite efficient calculations and keep the shapes still rather simple
224 times=[];
225 for ii = 1:length(grads)
226 g=grads{ii};
227 if is_trap(ii)
228 times = [times; cumsum([g.delay; g.riseTime; g.flatTime; g.fallTime])];
229 else
230 times = [times; g.delay+g.tt];
231 end
232 end
233 times=unique(times); % unique() also sorts the array
234 %times=unique(round(times/system.gradRasterTime)*system.gradRasterTime); % rounding to raster would be too crude here
235 dt=times(2:end)-times(1:end-1);
236 ieps=find(dt<eps);
237 if ~isempty(ieps)
238 dtx=[times(1); dt];
239 dtx(ieps)=dtx(ieps)+dtx(ieps+1); % this assumes that no more than two too similar values can occur
240 dtx(ieps+1)=[];
241 times=cumsum(dtx);
242 end
243 amplitudes=zeros(size(times));
244 for ii = 1:length(grads)
245 g=grads{ii};
246 if strcmp(g.type,'trap')
247 if g.flatTime>0 % trapezoid or triangle
248 g.tt=cumsum([0; g.riseTime; g.flatTime; g.fallTime]);
249 g.waveform=[0; g.amplitude; g.amplitude; 0];
250 else
251 g.tt=cumsum([0; g.riseTime; g.fallTime]);
252 g.waveform=[0; g.amplitude; 0];
253 end
254 end
255 tt=g.delay+g.tt;
256 % fix rounding for the first and last time points
257 [tmin, imin]=min(abs(tt(1)-times));
258 if tmin<eps
259 tt(1)=times(imin);
260 end
261 [tmin, imin]=min(abs(tt(end)-times));
262 if tmin<eps
263 tt(end)=times(imin);
264 end
265 % give up the "ownership" of the first point of the shape if it starts at a non-zero value
266 if abs(g.waveform(1))>eps && tt(1)>eps
267 tt(1)=tt(1)+eps;
268 end
269 amplitudes=amplitudes+interp1(tt,g.waveform,times,'linear',0);
270 end
271 grad=mr.makeExtendedTrapezoid(channel,'amplitudes',amplitudes,'times',times,'system',system);
272 return;
273end
275% OK, here we convert everything to a regularly-sampled waveform
276waveforms = {};
277max_length = 0;
278some_osa=any(is_osa);
279if some_osa
280 target_raster=system.gradRasterTime/2;
281else
282 target_raster=system.gradRasterTime;
283end
284for ii = 1:length(grads)
285 g = grads{ii};
286 if ~is_trap(ii)
287 if is_arb(ii)||is_osa(ii)
288 if some_osa && is_arb(ii)
289 % interpolate missing samples
290 waveforms{ii} = (g.waveform(floor(1:0.5:end))+g.waveform(ceil(1:0.5:end)))*0.5;
291 else
292 waveforms{ii} = g.waveform;
293 end
294 else
295 waveforms{ii} = mr.pts2waveform(g.tt, g.waveform, target_raster);
296 end
297 else
298 if (g.flatTime>0) % triangle or trapezoid
299 times = [g.delay - common_delay; ...
300 g.delay - common_delay + g.riseTime; ...
301 g.delay - common_delay + g.riseTime + g.flatTime; ...
302 g.delay - common_delay + g.riseTime + g.flatTime + g.fallTime];
303 amplitudes = [0; g.amplitude; g.amplitude; 0];
304 else
305 times = [g.delay - common_delay; ...
306 g.delay - common_delay + g.riseTime; ...
307 g.delay - common_delay + g.riseTime + g.fallTime];
308 amplitudes = [0; g.amplitude; 0];
309 end
310 waveforms{ii} = mr.pts2waveform(times, amplitudes, target_raster);
311 end
312 if size(waveforms{ii},1)==1
313 waveforms{ii}=waveforms{ii}';
314 end
315 %warning('addGradient(): potentially incorrect handling of delays... TODO: fixme!');
316 if g.delay - common_delay > 0
317 %warning('addGradient(): zerofilling the shape, running unchecked code...');
318 t_delay = (0:target_raster:g.delay-common_delay-target_raster).';
319 waveforms{ii} = [t_delay*0; waveforms{ii}];
320 end
321 max_length = max(max_length, length(waveforms{ii}));
322end
324w = zeros(max_length,1);
325for ii = 1:length(grads)
326 % % SK: Matlab is so ridiculously cumbersome...
327 % wt = zeros(max_length, 1);
328 % wt(1:length(waveforms{ii})) = waveforms{ii};
329 % w = w + wt;
330 % MZ: it is cumbersome indeed, but not so...
331 w(1:length(waveforms{ii})) = w(1:length(waveforms{ii})) + waveforms{ii};
332end
334grad = mr.makeArbitraryGrad(channel, w, system, ...
335 'maxSlew', maxSlew,...
336 'maxGrad', maxGrad,...
337 'delay', common_delay,...
338 'oversampling',some_osa,...
339 'first',sum(firsts(delays==common_delay)),...
340 'last',sum(lasts(durs==total_duration)));
341% first is defined by the sum of firsts with the minimal delay (common_delay)
342% last is defined by the sum of lasts with the maximum duration (total_duration)
343%grad.first=sum(firsts(delays==common_delay));
344%grad.last=sum(lasts(durs==total_duration));
346end