1function [grad, times, amplitudes] = makeExtendedTrapezoidArea(channel, grad_start, grad_end, area, sys)
2% Make the shortest possible extended trapezoid for a given area and edge values
3% This version is the one with the fixed flat top and was derived from the
4% corresponding PyPulseq version by Mehmet Emin Öztürk. This implementation
5% is both faster and more accurate than the previous one and it runs in Octave.
6% Main methodology in explained in Python version
7% Some variable names might also be different
9 if nargin < 5 || isempty(sys)
10 sys = default_opts(); % Define your default system
11 end
13 max_slew = sys.maxSlew * 0.99;
14 max_grad = sys.maxGrad * 0.99;
15 raster_time = sys.gradRasterTime;
17 min_duration = max(round(calc_ramp_time(grad_end, grad_start, max_slew, raster_time) / raster_time), 2);
19 % Estimate upper bound duration
20 max_duration = max( ...
21 [round(calc_ramp_time(0, grad_start, max_slew, raster_time) / raster_time), ...
22 round(calc_ramp_time(0, grad_end, max_slew, raster_time) / raster_time), ...
23 min_duration]);
25 % Try to find a solution linearly
26 solution = [];
27 for duration = min_duration:max_duration
28 solution = find_solution(duration, area, grad_start, grad_end, max_slew, max_grad, raster_time);
29 if ~isempty(solution)
30 break;
31 end
32 end
34 % Binary search if linear search failed
35 if isempty(solution)
36 duration = max_duration;
37 while isempty(solution)
38 duration = duration * 2;
39 solution = find_solution(duration, area, grad_start, grad_end, max_slew, max_grad, raster_time);
40 end
42 solution = binary_search(@(d) find_solution(d, area, grad_start, grad_end, max_slew, max_grad, raster_time), ...
43 floor(duration/2), duration);
44 end
46 time_ramp_up = solution(1) * raster_time;
47 flat_time = solution(2) * raster_time;
48 time_ramp_down = solution(3) * raster_time;
49 grad_amp = solution(4);
51 % Generate final time vector and amplitudes
52 if flat_time > 0
53 times = cumsum([0, time_ramp_up, flat_time, time_ramp_down]);
54 amplitudes = [grad_start, grad_amp, grad_amp, grad_end];
55 else
56 times = cumsum([0, time_ramp_up, time_ramp_down]);
57 amplitudes = [grad_start, grad_amp, grad_end];
58 end
60 grad=mr.makeExtendedTrapezoid(channel,'system',sys,'times',times, 'amplitudes', amplitudes);
61 if abs(grad.area - area) >= 1e-3
62 error('Could not find a solution for area=%.6f.', area);
63 end
64end
67function time = to_raster(time_val, raster_time)
68 time = ceil(time_val / raster_time) * raster_time;
69end
71function t = calc_ramp_time(g1, g2, max_slew, raster_time)
72 t = to_raster(abs(g1 - g2) / max_slew, raster_time);
73end
75function sol = binary_search(fun, low, high)
76 while low < high - 1
77 mid = floor((low + high) / 2);
78 if ~isempty(fun(mid))
79 high = mid;
80 else
81 low = mid;
82 end
83 end
84 sol = fun(high);
85end
88function sol = find_solution(duration, area, grad_start, grad_end, max_slew, max_grad, raster_time)
89 sign_area = sign(area);
90 if sign_area == 0
91 % zero-area request (e.g. linking two trajectory segments without
92 % changing the k-space position): sign(0) would zero out the search
93 % direction and produce Inf ranges in the flat=0 estimates below, so
94 % search on the side opposite to the edge gradients instead
95 sign_area = -sign(grad_start + grad_end);
96 if sign_area == 0
97 sign_area = 1;
98 end
99 end
100 grad_amp = sign_area * max_grad;
102 % Convert to raster steps
103 ru_min = abs(grad_amp - grad_start) / max_slew / raster_time;
104 rd_min = abs(grad_amp - grad_end) / max_slew / raster_time;
105 flat_time = max(duration - ru_min - rd_min, 0);
107 % Check if feasible
108 approx_area = ru_min * (grad_amp + grad_start) + ...
109 rd_min * (grad_amp + grad_end) + ...
110 2 * flat_time * grad_amp;
112 if abs(2 * area / raster_time) > abs(approx_area)
113 sol = [];
114 return;
115 end
117 % Easy early solution: max_grad
118 ru = (duration * max_slew * raster_time + sign_area * (grad_end - grad_start)) / (2 * max_slew * raster_time);
119 if sign_area * grad_start + ru * max_slew * raster_time > max_grad + 1e-5
120 ru_steps = round(abs(grad_start - sign_area * max_grad) / max_slew / raster_time);
121 rd_steps = round(abs(grad_end - sign_area * max_grad) / max_slew / raster_time);
122 flat_steps = duration - ru_steps - rd_steps;
123 if flat_steps > 0
124 grad_amp = -(ru_steps * raster_time * grad_start + ...
125 rd_steps * raster_time * grad_end - 2 * area) / ...
126 ((ru_steps + 2 * flat_steps + rd_steps) * raster_time);
127 amps = [grad_start, grad_amp, grad_amp, grad_end];
128 t = cumsum([0, ru_steps, flat_steps, rd_steps]) * raster_time;
129 slew = diff(amps) ./ diff(t);
130 if max(abs(slew)) < max_slew + 1e-5 && max(abs(amps)) < max_grad
131 sol = [ru_steps, flat_steps, rd_steps, grad_amp];
132 return;
133 end
134 end
135 end
137 % Conservative downscaling
138 while abs(2 * area / raster_time) < abs(approx_area)
139 grad_amp = grad_amp / 2;
140 if abs(grad_amp) < abs(max_grad) / 10
141 ru_min = 0; rd_min = 0;
142 flat_time = max(duration - ru_min - rd_min, 0);
143 break;
144 end
145 ru_min = abs(grad_amp - grad_start) / max_slew / raster_time;
146 rd_min = abs(grad_amp - grad_end) / max_slew / raster_time;
147 flat_time = max(duration - ru_min - rd_min, 0);
148 approx_area = ru_min * (grad_amp + grad_start) + ...
149 rd_min * (grad_amp + grad_end) + ...
150 2 * flat_time * grad_amp;
151 end
153 % Convert to integer steps
154 ru_min = floor(ru_min);
155 rd_min = floor(rd_min);
156 ru_limit = ceil(abs(sign_area * max_grad - grad_start) / max_slew / raster_time) + 1;
157 rd_limit = ceil(abs(sign_area * max_grad - grad_end) / max_slew / raster_time) + 1;
159 % All combinations
160 [RU, RD] = meshgrid(ru_min:ru_limit, rd_min:rd_limit);
161 RU = RU(:);
162 RD = RD(:);
163 valid_mask = RD < (duration - RU);
164 RU = RU(valid_mask);
165 RD = RD(valid_mask);
167 % Flat = 0 case
168 num = (2 * area - duration * (grad_end + grad_start) * raster_time);
169 denom_min = (grad_start - grad_end + sign_area * duration * max_slew * raster_time);
170 denom_max = (-grad_start + grad_end + sign_area * duration * max_slew * raster_time);
171 ru_flat0_min = round(num / denom_min / raster_time);
172 ru_flat0_max = duration - round(num / denom_max / raster_time);
174 RU = [RU; (ru_flat0_min:ru_flat0_max)'];
175 RD = [RD; (duration - (ru_flat0_min:ru_flat0_max))'];
177 % Filter invalid ones
178 flat = duration - RU - RD;
179 valid = flat >= 0 & RU > 0 & RD > 0;
180 RU = RU(valid); RD = RD(valid); flat = flat(valid);
182 % Calculate amp
183 grad_amp = -(RU * raster_time * grad_start + RD * raster_time * grad_end - 2 * area) ./ ...
184 ((RU + 2 * flat + RD) * raster_time);
186 % Slew
187 slew1 = abs(grad_start - grad_amp) ./ (RU * raster_time);
188 slew2 = abs(grad_end - grad_amp) ./ (RD * raster_time);
190 % Valid gradient/slew combinations
191 valid = abs(grad_amp) <= max_grad + 1e-5 & slew1 <= max_slew + 1e-5 & slew2 <= max_slew + 1e-5;
192 if ~any(valid)
193 sol = [];
194 return;
195 end
197 % Pick lowest slew
198 ind = find(valid);
199 [~, min_idx] = min(slew1(ind) + slew2(ind));
200 best = ind(min_idx);
202 sol = [RU(best), flat(best), RD(best), grad_amp(best)];
203end