78d04f1seqlab: write and view pulseq MRI sequences in the browserJeremy Magland 1function [kout, success] = calcRamp(k0,kend,varargin)
2%
3% the aim of joink is to join the points k0 and kend in three-dimensional
4% k-space in minimal time, observing the gradient and slew limits, and the
5% gradient strength G0 before k0(:,2) and Gend after kend(:,1)
6%
7% In the context of a fixed gradient dwell time this is a discrete problem
8% with an a priori unknown number of discretization steps. Therefore joink
9% tries out the optimization with 0 steps, then 1 step, and so on, until
10% all conditions can be fulfilled, thus yielding a short connection
11%
12% N.B. The connection found this way is not necessarily always the shortest
13% (there are some counterexamples) but still quite short. Improvements
14% possible.
15%
16% Usage: [kout success] = joink(k0,kend,MaxGrad,MaxSlew,GradDwell,MaxPoints)
17%
18% [kout] connecting k-space points without k0 and kend, size = [3,Nt],
19% where Nt = number of steps between k0 and kend.
20% k-space units: 1/m
21%
22% [success] a flag indicating if a solution was found with up to
23% MaxPoints k-space points: success (1), no solution (0)
24%
25% [k0] Two preceding points in k-space, size = [3,2]. From these
26% points, the starting gradient will be calculated.
27%
28% [kend] Two following points in k-space, size = [3,2]. From these
29% points, the target gradient will be calculated.
30%
31% [MaxGrad] maximum total vector gradient strength, size = [1,1]
32% or
33% [MaxGrad] maximum gradient strength per coordinate, size = [3,1]
34%
35% [MaxSlew] maximum total vector slew rate, size = [1,1],
36% or
37% [MaxSlew] maximum slew rate per coordinate, size = [3,1]
38% all slew units: T/(m*s)
39%
40% [GradDwell] time between two k-space points, size = [1,1], unit: s
41%
42% [MaxPoints] maximum number of k-space points to be used in connecting k0
43% with kend. Keep at a reasonable order of magnitude!
44%
46persistent parser
47if isempty(parser)
48 parser = mr.aux.InputParserCompat;
49 parser.FunctionName = 'calcRamp';
50 parser.addRequired('k0',@isnumeric);
51 parser.addRequired('kend',@isnumeric);
52 parser.addOptional('system',[],@isstruct);
53 parser.addParamValue('MaxPoints',500,@isnumeric);
54 parser.addParamValue('maxGrad',0,@isnumeric);
55 parser.addParamValue('maxSlew',0,@isnumeric);
56 parser.addParamValue('gradOversampling',false,@islogical);
57end
58parse(parser,k0,kend,varargin{:});
59opt = parser.Results;
61if isempty(opt.system)
62 system=mr.opts();
63else
64 system=opt.system;
65end
67maxSlew=system.maxSlew;
68maxGrad=system.maxGrad;
69if opt.maxGrad>0
70 maxGrad=opt.maxGrad;
71end
72if opt.maxSlew>0
73 maxSlew=opt.maxSlew;
74end
75if opt.gradOversampling
76 GradRaster = system.gradRasterTime/2;
77else
78 GradRaster = system.gradRasterTime;
79end
80MaxPoints = opt.MaxPoints;
82if ~mr.aux.isOctave()
83 SaveRecLimit = get(0,'RecursionLimit');
84 set(0,'RecursionLimit',MaxPoints+10);
85end
87% Determine whether we are in componentwise limited mode or in total vector
88% limited mode.
89% mode = 0 --> total vector limited
90% mode = 1 --> componentwise limited
92if isequal(size(maxGrad),[1 1]) && isequal(size(maxSlew),[1 1])
93 mode = 0;
94elseif isequal(size(maxGrad),[3 1]) && isequal(size(maxSlew),[3 1])
95 mode = 1;
96else
97 error('Input value MaxGrad or MaxSlew in invalid format.');
98end
100G0 = (k0 (:,2)-k0 (:,1))/GradRaster;
101Gend = (kend(:,2)-kend(:,1))/GradRaster;
102k0 = k0 (:,2);
103kend = kend(:,1);
105success = 0;
106kout = zeros(3,0);
107UsePoints = 0; % first try: connecting directly
109while (success == 0) && (UsePoints <= MaxPoints)
110 if mode == 0
111 if (norm(G0)>maxGrad) || (norm(Gend)>maxGrad)
112 break;
113 end;
114 kout = joinleft0(k0,kend,G0,Gend,UsePoints);
115 else
116 if (abs(G0)>abs(maxGrad)) || (abs(Gend)>abs(maxGrad))
117 break;
118 end;
119 kout = joinleft1(k0,kend,G0,Gend,UsePoints);
120 end
121 UsePoints = UsePoints + 1;
122end
124if ~mr.aux.isOctave()
125 set(0,'RecursionLimit',SaveRecLimit); % set previous recursion limit
126end
128% -------------------------------------------------------------------------
129function ok = InsideLimits(Grad,Slew)
130%
131% check if both gradient and slew rates are inside the respective limits
132%
133if mode == 0
134 Grad2 = sum(Grad.^2,1); % gradient vector norm squared
135 Slew2 = sum(Slew.^2,1); % slew vector norm squared
136 ok = (max(Grad2) <= maxGrad^2) && (max(Slew2) <= maxSlew^2);
137else
138 ok = (sum(max(abs(Grad),[],2) <= maxGrad) == 3) && ...
139 (sum(max(abs(Slew),[],2) <= maxSlew) == 3);
140end
142end % function InsideLimits
144% -------------------------------------------------------------------------
145function koutleft = joinleft0(k0,kend,G0,Gend,UsePoints)
146%
147% Add one k-space point close to k0. Gradient and slew limits apply in
148% total vector limited mode.
149%
150% Rationale:
151%
152% 0. If UsePoints == 0 the recursion stops. If k0 and kend can be joined
153% in one GradDwell time, return success, else return "no success".
154%
155% 1. Calculate optimal k-space point kopt that would lie on a straight
156% line of N=UsePoints evenly spaced points to kend. If this kopt can be
157% reached within gradient and slew limts, kopt is the solution of this
158% function call.
159%
160% 2. If kopt cannot be reached, calculate the gradient limited point kgl
161% closest to kopt. If this point can be reached in one GradDwell time
162% without violating the slew limit, kgl is the solution of this
163% function call.
164%
165% 3. If not kgl is not inside the slew limit, the slew limited point
166% closest to kopt, ksl, is calculated. If ksl is inside the gradient
167% limit, ksl is the solution.
168%
169% 4. If neither kgl nor ksl are possible find the point kglsl closest to
170% kopt that satisfies both limits at the same time. See
171% illustration.fig / illustration.png
172%
173% 5. Call joinright0 to obtain the other points starting with a point
174% next to kend.
175%
176if UsePoints == 0
177 G = [G0 (kend-k0)/GradRaster Gend];
178 S = (G(:,2:end)-G(:,1:end-1))/GradRaster;
180 koutleft = zeros(3,0); % no additional k-space point
181 success = InsideLimits(G,S);
183 return;
184end
186dk = (kend-k0)/(UsePoints+1);
187kopt = k0+dk; % this would be on the direct
188Gopt = (kopt-k0)/GradRaster; % line
189Sopt = (Gopt-G0)/GradRaster;
191okGopt = (sum(Gopt.^2,1) <= maxGrad^2);
192okSopt = (sum(Sopt.^2,1) <= maxSlew^2);
194if okGopt && okSopt
195 kLeft = kopt;
196else
197 a = maxGrad*GradRaster;
198 b = maxSlew*GradRaster^2;
200 dkprol = G0*GradRaster; % prolonged point with no change
201 dkconn = dk-dkprol; % in gradient
202 % slew limited closest to kopt
203 ksl = k0 + dkprol + dkconn/norm(dkconn)*b;
204 Gsl = (ksl-k0)/GradRaster;
205 okGsl = (sum(Gsl.^2,1) <= maxGrad^2);
207 kgl = k0 + dk/norm(dk)*a; % gradient limited closest to
208 Ggl = (kgl-k0)/GradRaster; % kopt
209 Sgl = (Ggl-G0)/GradRaster;
210 okSgl = (sum(Sgl.^2,1) <= maxSlew^2);
212 if okGsl
213 kLeft = ksl;
214 elseif okSgl
215 kLeft = kgl;
216 else
217 c = norm(dkprol);
218 c1 = (a^2-b^2+c^2)/(2*c); % if e.g. |Gend|<Gmax, then
219 h = sqrt(a^2-c1^2); % a^2-c1^2 is always positive
220 kglsl = k0 + c1*dkprol/norm(dkprol);
221 projondkprol = (kgl*dkprol.') * dkprol/norm(dkprol);
222 hdirection = kgl - projondkprol;
223 kglsl = kglsl + h*hdirection/norm(hdirection);
224 kLeft = kglsl;
225 end
226end
228k = joinright0(kLeft,kend,(kLeft-k0)/GradRaster,Gend,UsePoints-1);
230koutleft = [kLeft k]; % pass along result
232end % function joinleft0
234% -------------------------------------------------------------------------
235function koutright = joinright0(k0,kend,G0,Gend,UsePoints)
236%
237% Add one k-space point close to kend. Gradient and slew limits apply in
238% total vector limited mode. Rationale see joinleft0.
239%
240if UsePoints == 0
241 G = [G0 (kend-k0)/GradRaster Gend];
242 S = (G(:,2:end)-G(:,1:end-1))/GradRaster;
244 koutright = zeros(3,0); % no additional k-space point
245 success = InsideLimits(G,S);
247 return;
248end
250dk = (k0-kend)/(UsePoints+1);
251kopt = kend+dk; % this would be on the direct
252Gopt = (kend-kopt)/GradRaster; % line
253Sopt = (Gend-Gopt)/GradRaster;
255okGopt = (sum(Gopt.^2,1) <= maxGrad^2);
256okSopt = (sum(Sopt.^2,1) <= maxSlew^2);
258if okGopt && okSopt
259 kRight = kopt;
260else
261 a = maxGrad*GradRaster;
262 b = maxSlew*GradRaster^2;
264 dkprol = -Gend*GradRaster; % prolonged point with no change
265 dkconn = dk-dkprol; % in gradient
266 % slew limited closest to kopt
267 ksl = kend + dkprol + dkconn/norm(dkconn)*b;
268 Gsl = (kend-ksl)/GradRaster;
269 okGsl = (sum(Gsl.^2,1) <= maxGrad^2);
271 kgl = kend + dk/norm(dk)*a; % gradient limited closest to
272 Ggl = (kend-kgl)/GradRaster; % kopt
273 Sgl = (Gend-Ggl)/GradRaster;
274 okSgl = (sum(Sgl.^2,1) <= maxSlew^2);
276 if okGsl
277 kRight = ksl;
278 elseif okSgl
279 kRight = kgl;
280 else
281 c = norm(dkprol);
282 c1 = (a^2-b^2+c^2)/(2*c); % if e.g. |Gend|<Gmax, then
283 h = sqrt(a^2-c1^2); % a^2-c1^2 is always positive
284 kglsl = kend + c1*dkprol/norm(dkprol);
285 projondkprol = (kgl*dkprol.') * dkprol/norm(dkprol);
286 hdirection = kgl - projondkprol;
287 kglsl = kglsl + h*hdirection/norm(hdirection);
288 kRight = kglsl;
289 end
290end
292k = joinleft0(k0,kRight,G0,(kend-kRight)/GradRaster,UsePoints-1);
294koutright = [k kRight]; % pass along result
296end % function joinright0
298% -------------------------------------------------------------------------
299function koutleft = joinleft1(k0,kend,G0,Gend,UsePoints)
300%
301% Add one k-space point close to k0. Gradient and slew limits apply in
302% componentwise limited mode. Rationale is the same as in joinleft0 but
303% it's much easier to find the point kglsl in step 4.
304%
305if UsePoints == 0
306 G = [G0 (kend-k0)/GradRaster Gend];
307 S = (G(:,2:end)-G(:,1:end-1))/GradRaster;
309 koutleft = zeros(3,0); % no additional k-space point
310 success = InsideLimits(G,S);
312 return;
313end
315kLeft = zeros(3,1);
317dk = (kend-k0)/(UsePoints+1);
318kopt = k0+dk; % this would be on the direct
319Gopt = (kopt-k0)/GradRaster; % line
320Sopt = (Gopt-G0)/GradRaster;
321okGopt = (abs(Gopt) <= maxGrad);
322okSopt = (abs(Sopt) <= maxSlew);
324dkprol = G0*GradRaster;
325dkconn = dk-dkprol;
326 % slew limited
327ksl = k0 + dkprol + sign(dkconn).*maxSlew*GradRaster^2;
328Gsl = (ksl-k0)/GradRaster;
329okGsl = (abs(Gsl) <= maxGrad);
330 % gradient limited
331kgl = k0 + sign(dk).*maxGrad*GradRaster;
332Ggl = (kgl-k0)/GradRaster;
333Sgl = (Ggl-G0)/GradRaster;
334okSgl = (abs(Sgl) <= maxSlew);
336for ii=1:3
337 if (okGopt(ii)==1) && (okSopt(ii)==1)
338 kLeft(ii) = kopt(ii);
339 elseif (okGsl(ii)==1)
340 kLeft(ii) = ksl(ii);
341 elseif (okSgl(ii)==1)
342 kLeft(ii) = kgl(ii);
343 else
344 display('Moment mal - hier dürfte ich niemals hinkommen!');
345 end
346end
348k = joinright1(kLeft,kend,(kLeft-k0)/GradRaster,Gend,UsePoints-1);
350koutleft = [kLeft k]; % pass along result
352end % function joinleft1
354% -------------------------------------------------------------------------
355function koutright = joinright1(k0,kend,G0,Gend,UsePoints)
356%
357% Add one k-space point close to kend. Gradient and slew limits apply in
358% componentwise limited mode. Rationale is the same as in joinright1
359%
360if UsePoints == 0
361 G = [G0 (kend-k0)/GradRaster Gend];
362 S = (G(:,2:end)-G(:,1:end-1))/GradRaster;
364 koutright = zeros(3,0); % no additional k-space point
365 success = InsideLimits(G,S);
367 return;
368end
370kRight = zeros(3,1);
372dk = (k0-kend)/(UsePoints+1);
373kopt = kend+dk; % this would be on the direct
374Gopt = (kend-kopt)/GradRaster; % line
375Sopt = (Gend-Gopt)/GradRaster;
376okGopt = (abs(Gopt) <= maxGrad);
377okSopt = (abs(Sopt) <= maxSlew);
379dkprol = -Gend*GradRaster;
380dkconn = dk-dkprol;
381 % slew limited
382ksl = kend + dkprol + sign(dkconn).*maxSlew*GradRaster^2;
383Gsl = (kend-ksl)/GradRaster;
384okGsl = (abs(Gsl) <= maxGrad);
385 % gradient limited
386kgl = kend + sign(dk).*maxGrad*GradRaster;
387Ggl = (kend-kgl)/GradRaster;
388Sgl = (Gend-Ggl)/GradRaster;
389okSgl = (abs(Sgl) <= maxSlew);
391for ii = 1:3
392 if (okGopt(ii)==1) && (okSopt(ii)==1)
393 kRight(ii) = kopt(ii);
394 elseif (okGsl(ii)==1)
395 kRight(ii) = ksl(ii);
396 elseif (okSgl(ii)==1)
397 kRight(ii) = kgl(ii);
398 else
399 error('Unknown error. Code should not execute');
400 end
401end
403k = joinleft1(k0,kRight,G0,(kend-kRight)/GradRaster,UsePoints-1);
405koutright = [k kRight]; % pass along result
407end % function joinright1
409end % main
410% =========================================================================