1function [grad, times, amplitudes] = makeHexagonGradientArea(channel, grad_start, grad_end, area, sys)
2% Make the shortest possible hexagonal gradien (creating an extended
3% trapezoid object) for a given area and edge values. In contrast to
4% mr.makeExtendedTrapezoidArea(), this function creates a generic
5% polynomial gradient object without a plato between the vortex2 and
6% vortex3. We expect this function to find better solutions by relaxing the
7% constraint of two vertices having the same amplitude, however, as seen in
8% testCase_12, this objective cannot be achieved in all cases by the current
9% implementation.
10% Parameters and methodology are explained in mr.makeExtendedTrapezoidArea().
11% This function was implemented by Mehmet Emin Öztürk during his visit to
12% Freiburg with some input from Maxim Zaitsev.
14 if nargin < 5 || isempty(sys)
15 sys = default_opts(); % Define your default system
16 end
18 % validate inputs the search below cannot handle (runaway loops or obscure downstream errors)
19 if abs(grad_start) > sys.maxGrad
20 error('grad_start amplitude violation (%.0f%%)', abs(grad_start)/sys.maxGrad*100);
21 end
22 if abs(grad_end) > sys.maxGrad
23 error('grad_end amplitude violation (%.0f%%)', abs(grad_end)/sys.maxGrad*100);
24 end
26 max_slew = sys.maxSlew * 0.99;
27 max_grad = sys.maxGrad * 0.99;
28 raster_time = sys.gradRasterTime;
30 min_duration = max(round(calc_ramp_time(grad_end, grad_start, max_slew, raster_time) / raster_time), 2);
32 % Estimate upper bound duration
33 max_duration = max( ...
34 [round(calc_ramp_time(0, grad_start, max_slew, raster_time) / raster_time), ...
35 round(calc_ramp_time(0, grad_end, max_slew, raster_time) / raster_time), ...
36 min_duration]);
38 % Try to find a solution linearly
39 times = [];
40 amplitudes = [];
41 for duration = min_duration:max_duration
42 [times, amplitudes] = find_solution(duration, area, grad_start, grad_end, max_slew, max_grad, raster_time);
43 if ~isempty(times)
44 break;
45 end
46 end
48 % Binary search if linear search failed
49 if isempty(times)
50 duration = max_duration;
51 while isempty(times)
52 duration = duration * 2;
53 [times, ~] = find_solution(duration, area, grad_start, grad_end, max_slew, max_grad, raster_time);
54 end
56 [times, amplitudes] = binary_search(@(d) find_solution(d, area, grad_start, grad_end, max_slew, max_grad, raster_time), ...
57 floor(duration/2), duration);
58 end
59 % drop zero-length segments (their end points always carry equal
60 % amplitudes, otherwise the slew-rate check would have rejected them)
61 keep = [true, diff(times) > 0];
62 times = times(keep);
63 amplitudes = amplitudes(keep);
64 grad=mr.makeExtendedTrapezoid(channel,'system',sys,'times',times, 'amplitudes', amplitudes);
65 if abs(grad.area - area) >= 1e-3
66 error('Could not find a solution for area=%.6f.', area);
67 end
68end
71function time = to_raster(time_val, raster_time)
72 time = ceil(time_val / raster_time) * raster_time;
73end
75function t = calc_ramp_time(g1, g2, max_slew, raster_time)
76 t = to_raster(abs(g1 - g2) / max_slew, raster_time);
77end
79function [times, amplitudes] = binary_search(fun, low, high)
80 while low < high - 1
81 mid = floor((low + high) / 2);
82 if ~isempty(fun(mid))
83 high = mid;
84 else
85 low = mid;
86 end
87 end
88 [times, amplitudes] = fun(high);
89end
92function [times, amplitudes] = find_solution(duration, area, grad_start, grad_end, max_slew, max_grad, raster_time)
93% Find extended trapezoid gradient waveform for given duration
95 sign_area = sign(area);
96 if sign_area == 0
97 % zero-area request (e.g. linking two trajectory segments without
98 % changing the k-space position): sign(0) would zero out the search
99 % direction and send the loops below into runaway iterations, so
100 % search on the side opposite to the edge gradients instead
101 sign_area = -sign(grad_start + grad_end);
102 if sign_area == 0
103 sign_area = 1;
104 end
105 end
106 grad_amp = sign_area * max_grad;
107 ramp_up_times = [];
108 ramp_down_times = [];
110 % Early estimation
111 ru_min = abs(grad_amp - grad_start) / max_slew / raster_time;
112 rd_min = abs(grad_amp - grad_end) / max_slew / raster_time;
113 flat_time = max(duration - ru_min - rd_min, 0);
115 area_check = ru_min * (grad_amp + grad_start) + rd_min * (grad_amp + grad_end) + 2 * flat_time * grad_amp;
117 if abs(2 * area / raster_time) > abs(area_check)
118 times = []; amplitudes = [];
119 return;
120 end
122 % Fast-case: max_grad ramp with valid flat
123 ru = (duration * max_slew * raster_time + sign_area * (grad_end - grad_start)) / (2 * max_slew * raster_time);
124 if sign_area * grad_start + ru * max_slew * raster_time > max_grad + 1e-5
125 ru_steps = round(abs(grad_start - sign_area * max_grad) / max_slew / raster_time);
126 rd_steps = round(abs(grad_end - sign_area * max_grad) / max_slew / raster_time);
127 flat_steps = duration - ru_steps - rd_steps;
128 if flat_steps > 0
129 grad_amp = -(ru_steps * raster_time * grad_start + rd_steps * raster_time * grad_end - 2 * area) / ...
130 ((ru_steps + 2 * flat_steps + rd_steps) * raster_time);
131 amps = [grad_start, grad_amp, grad_amp, grad_end];
132 t = cumsum([0, ru_steps, flat_steps, rd_steps]) * raster_time;
133 slew = diff(amps) ./ diff(t);
134 if max(abs(slew)) < max_slew + 1e-5 && max(abs(amps)) < max_grad
135 times = t;
136 amplitudes = amps;
137 return;
138 end
139 end
140 end
142 % Gradually reduce grad_amp if area too large
143 while abs(2 * area / raster_time) < abs(area_check)
144 grad_amp = grad_amp / 2;
145 if abs(grad_amp) < abs(max_grad) / 10
146 ru_min = 0; rd_min = 0;
147 flat_time = max(duration - ru_min - rd_min, 0);
148 break;
149 end
150 ru_min = abs(grad_amp - grad_start) / max_slew / raster_time;
151 rd_min = abs(grad_amp - grad_end) / max_slew / raster_time;
152 flat_time = max(duration - ru_min - rd_min, 0);
153 area_check = ru_min * (grad_amp + grad_start) + rd_min * (grad_amp + grad_end) + 2 * flat_time * grad_amp;
154 end
156 % Discrete timing limits
157 ru_min = floor(ru_min);
158 rd_min = floor(rd_min);
159 ru_limit = ceil(abs(sign_area * max_grad - grad_start) / max_slew / raster_time);
160 rd_limit = ceil(abs(sign_area * max_grad - grad_end) / max_slew / raster_time);
162 flat_time = duration - min(rd_min, rd_limit) - min(ru_min, ru_limit);
163 flat_time_min = duration - rd_limit - ru_limit;
165 min_dif_area = area;
166 i = -1;
168 while flat_time > max(flat_time_min, -1)
169 i = i + 1;
170 ru_max = ru_min + i;
171 rd_max = rd_min + i;
172 flat_time = duration - min(rd_min + i, rd_limit) - min(ru_min + i, ru_limit);
174 if flat_time <= 0
175 % Fallback: flat_time = 0 or 1
176 ru_0min = (2 * area - duration * (grad_end + grad_start) * raster_time) / ...
177 (grad_start - grad_end + sign_area * duration * max_slew * raster_time) / raster_time;
178 ru_0max = duration - (2 * area - duration * (grad_end + grad_start) * raster_time) / ...
179 (-grad_start + grad_end + sign_area * duration * max_slew * raster_time) / raster_time;
181 ru_0min = floor(ru_0min); ru_0max = ceil(ru_0max);
182 for ru_try = ru_0min:ru_0max
183 ramp_up_times = [ramp_up_times, ru_try, ru_try];
184 ramp_down_times = [ramp_down_times, duration - ru_try - 1, duration - ru_try];
185 end
186 break;
187 end
189 % Gradients at corners
190 grad_p0 = sign_area * max_slew * ru_max * raster_time + grad_start;
191 grad_p1 = sign_area * max_slew * rd_max * raster_time + grad_end;
193 if abs(grad_p0) >= max_grad
194 limit_option1 = abs(sign_area * max_slew * (ru_limit - 1) * raster_time + grad_start);
195 limit_option2 = abs(((ru_limit + flat_time) * sign_area * max_grad + grad_start - grad_p1) / (ru_limit + flat_time));
196 if limit_option1 > limit_option2
197 grad_p0 = sign_area * max_slew * (ru_limit - 1) * raster_time + grad_start;
198 ru_max = floor(abs(sign_area * max_grad - grad_start) / max_slew / raster_time);
199 ru_limit = ru_max;
200 else
201 grad_p0 = sign_area * max_grad;
202 ru_max = ceil(abs(sign_area * max_grad - grad_start) / max_slew / raster_time);
203 ru_limit = ru_max;
204 end
205 flat_time = duration - rd_max - ru_max;
206 end
208 if abs(grad_p1) >= max_grad
209 limit_option1 = abs(sign_area * max_slew * (rd_limit - 1) * raster_time + grad_end);
210 limit_option2 = abs(((rd_limit + flat_time) * sign_area * max_grad + grad_end - grad_p0) / (rd_limit + flat_time));
211 if limit_option1 > limit_option2
212 grad_p1 = sign_area * max_slew * (rd_limit - 1) * raster_time + grad_end;
213 rd_max = floor(abs(sign_area * max_grad - grad_end) / max_slew / raster_time);
214 rd_limit = rd_max;
215 else
216 grad_p1 = sign_area * max_grad;
217 rd_max = ceil(abs(sign_area * max_grad - grad_end) / max_slew / raster_time);
218 rd_limit = rd_max;
219 end
220 flat_time = duration - rd_max - ru_max;
221 end
223 % Slope too high from p0 to p1
224 if abs(grad_p0 - grad_p1) / (flat_time * raster_time) > max_slew
225 if abs(grad_p0) < abs(grad_p1)
226 grad_p1 = grad_p0 + flat_time * raster_time * max_slew * sign_area;
227 rd_max = ceil(abs(grad_end - grad_p1) / max_slew / raster_time);
228 else
229 grad_p0 = grad_p1 + flat_time * raster_time * max_slew * sign_area;
230 ru_max = ceil(abs(grad_start - grad_p0) / max_slew / raster_time);
231 end
232 flat_time = duration - rd_max - ru_max;
233 end
235 % Compute area
236 area_current = raster_time/2 * ...
237 (ru_max * (grad_p0 + grad_start) + ...
238 flat_time * (grad_p0 + grad_p1) + ...
239 rd_max * (grad_p1 + grad_end));
241 if sign(min_dif_area) ~= sign(area - area_current)
242 t = round(cumsum([0, ru_max, flat_time, rd_max]) * raster_time, 5);
243 if abs(grad_p0) < abs(grad_p1)
244 Gtest = grad_p1;
245 corner_grad = -(grad_start * ru_max + grad_end * rd_max + Gtest * (flat_time + rd_max) - 2 * area / raster_time) / (flat_time + ru_max);
246 amps = [grad_start, corner_grad, Gtest, grad_end];
247 else
248 Gtest = grad_p0;
249 corner_grad = -(grad_start * ru_max + grad_end * rd_max + Gtest * (flat_time + ru_max) - 2 * area / raster_time) / (flat_time + rd_max);
250 amps = [grad_start, Gtest, corner_grad, grad_end];
251 end
253 if any(round(diff(t)/raster_time) < 1)
254 continue;
255 end
257 slew = diff(amps) ./ diff(t);
258 if max(abs(slew)) <= max_slew + 1e-5
259 times = t;
260 amplitudes = amps;
261 return;
262 end
264 % fallback search if slope still too large
265 for ru_try = ru_max-1 : duration - rd_max
266 for rd_try = rd_max-1 : duration - ru_try
267 flat = duration - ru_try - rd_try;
268 grad_amp = -(ru_try * raster_time * grad_start + rd_try * raster_time * grad_end - 2 * area) / ...
269 ((ru_try + 2 * flat + rd_try) * raster_time);
270 amps = [grad_start, grad_amp, grad_amp, grad_end];
271 t = cumsum([0, ru_try, flat, rd_try]) * raster_time;
272 slew = diff(amps) ./ diff(t);
273 if max(abs(slew)) < max_slew + 1e-5 && max(abs(amps)) < max_grad
274 times = t;
275 amplitudes = amps;
276 return;
277 end
278 end
279 end
280 end
282 if abs(area_current - area) < abs(min_dif_area)
283 min_dif_area = area - area_current;
284 end
285 end
287 % Fallback to triangle search
288 ru_vec = ramp_up_times(:);
289 rd_vec = ramp_down_times(:);
290 valid = ru_vec .* rd_vec > 0;
291 ru_vec = ru_vec(valid); rd_vec = rd_vec(valid);
292 flat = duration - ru_vec - rd_vec;
293 valid = flat >= 0;
294 ru_vec = ru_vec(valid); rd_vec = rd_vec(valid); flat = flat(valid);
296 grad_amp = -(ru_vec * raster_time * grad_start + rd_vec * raster_time * grad_end - 2 * area) ./ ...
297 ((ru_vec + 2 * flat + rd_vec) * raster_time);
299 slew1 = abs(grad_start - grad_amp) ./ (ru_vec * raster_time);
300 slew2 = abs(grad_end - grad_amp) ./ (rd_vec * raster_time);
302 valid = abs(grad_amp) <= max_grad + 1e-5 & slew1 <= max_slew + 1e-5 & slew2 <= max_slew + 1e-5;
303 idx = find(valid, 1);
304 if isempty(idx)
305 times = []; amplitudes = [];
306 return;
307 end
309 t = cumsum([0, ru_vec(idx), flat(idx), rd_vec(idx)]) * raster_time;
310 amps = [grad_start, grad_amp(idx), grad_amp(idx), grad_end];
312 times = t;
313 amplitudes = amps;
314end