JEMRIS 2.9.2
open-source MRI simulations
Loading...
Searching...
No Matches
mtg_functions.h
1#ifndef MTG_FUNCTIONS
2#define MTG_FUNCTIONS 1
3
4/*This file contains all the functions used to calculated the time optimal gradient waveforms.
5
6minTimeGradientRIV - Computes the rotationally invariant solution
7 RungeKutte_riv - Used to solve the ODE using RK4
8 beta - calculates sqrt (gamma^2 * smax^2 - k^2 * st^4) in the ODE. Used in RungeKutte_riv
9minTimeGradientRV - Computes the rotationally variant solution
10 RungeKutte_rv - Used to solve the ODE using RK4
11 sdotdot - calculates the maximum possible value for d^2s/dt^t, used in RugeKutte_rv */
12
13#include "spline.h"
14#include "float.h"
15#include "math.h"
16
17inline double beta(double k, double st, double smax, double gamma = 4.25756) {
18 /* calculates sqrt (gamma^2 * smax^2 - k^2 * st^4) used in RK4 method for rotationally invariant ODE solver */
19 return sqrt(sqrt((gamma*gamma*smax*smax - k*k*st*st*st*st)*(gamma*gamma*smax*smax - k*k*st*st*st*st)));
20}
21
22inline double RungeKutte_riv(double ds, double st, double k[], double smax, double gamma = 4.25756) {
23 /* Solves ODE for rotationally invariant solution using Runge-Kutte method*/
24 double k1 = ds * (1/st) * beta(k[0], st, smax, gamma);
25 double k2 = ds * 1 / (st + k1/2) * beta(k[1], st + k1/2, smax, gamma);
26 double k3 = ds * 1 / (st + k2/2) * beta(k[1], st + k2/2, smax, gamma);
27 double k4 = ds * 1 / (st + k3/2) * beta(k[2], st + k3/2, smax, gamma);
28 return (k1/6 + k2/3 + k3/3 + k4/6);
29}
30
31inline void minTimeGradientRIV(double* x, double* y, double* z, int Lp, double g0, double gfin, double gmax, double smax,
32 double dt, double*& gx, double*& gy, double*& gz,
33 int &l_t, double ds = -1, double gamma = 4.25756) {
34
35 /*Finds the time optimal gradient waveforms for the rotationally invariant constraints case.
36
37 x - The input curve in x
38 y - , y
39 z - , and z dir.
40 Lp - row dimension of x,y,z
41 g0 - Initial gradient amplitude.
42 gfin - Gradient value at the end of the trajectory.
43 If given value is not possible
44 the result would be the largest possible amplitude.
45 gmax - Maximum gradient [G/cm] (4 default)
46 smax - Maximum slew [G/cm/ms] (15 default)
47 dt - Sampling time intervale [ms] (4e-3 default)
48 gx, gy, gz - pointers to gradient waveforms to be returned
49 kx, ky, kz - pointers to k-space trajectory after interpolation
50 l_t - Dimension of interpolated k-space trajectory (kx, ky, kz) needed for creating mex return arrays.
51 */
52
53 int i = 0;
54
55 /* iflag used in spline method to signal error */
56 int iflag;
57 int last;
58
59 double *p = new double[Lp];
60
61 /* Representing the curve with parameter p */
62 for (i = 0; i < Lp; ++i)
63 p[i] = i;
64
65 /* Interpolation of curve for gradient accuracy, using cubic spline interpolation */
66
67 /* arrays used by spline function to store coefficients. */
68 double *c1x = new double[Lp];
69 double *c2x = new double[Lp];
70 double *c3x = new double[Lp];
71 double *c1y = new double[Lp];
72 double *c2y = new double[Lp];
73 double *c3y = new double[Lp];
74 double *c1z = new double[Lp];
75 double *c2z = new double[Lp];
76 double *c3z = new double[Lp];
77
78 spline(Lp, 0, 0, 1, 1, p, x, c1x, c2x, c3x, &iflag);
79 spline(Lp, 0, 0, 1, 1, p, y, c1y, c2y, c3y, &iflag);
80 spline(Lp, 0, 0, 1, 1, p, z, c1z, c2z, c3z, &iflag);
81
82 double dp = 0.1;
83 int num_evals = (int) floor((Lp-1) / dp)+ 1;
84
85 double toeval = 0;
86
87 /* converting to arc-length parameterization from p, using trapezoidal integration */
88 double *s_of_p = new double[num_evals];
89 s_of_p[0] = 0;
90
91 double *sop_num = new double[num_evals];
92
93 double sofar = 0;
94 double Cp_abs_pre = 0.;
95 for (i = 0; i < num_evals; ++i) {
96 toeval = (double) i * dp;
97 double Cpx = deriv(Lp, toeval, p, c1x, c2x, c3x, &last);
98 double Cpy = deriv(Lp, toeval, p, c1y, c2y, c3y, &last);
99 double Cpz = deriv(Lp, toeval, p, c1z, c2z, c3z, &last);
100 /* interpolated curve in p-parameterization */
101 double Cp_abs = sqrt(Cpx*Cpx + Cpy*Cpy + Cpz*Cpz);
102 sofar += (Cp_abs + Cp_abs_pre)/2;
103 s_of_p[i] = dp * sofar;
104 Cp_abs_pre = Cp_abs;
105
106 sop_num[i] = i * dp;
107 }
108
109 /* length of the curve */
110 double L = s_of_p[num_evals-1];
111
112 /* decide ds and compute st for the first point */
113 double stt0 = gamma*smax; /* always assumes first point is max slew */
114 double st0 = (stt0*dt)/2; /* start at half the gradient for accuracy close to g=0 */
115 double s0 = st0*dt;
116
117 if (ds < 0) { /* if a ds value was not specified */
118// ds = s0/1.5; /* smaller step size for numerical accuracy */
119 ds = fabs(ds) * s0/1.5; /* pehses: to make it easier to scale ds */
120 }
121
122 int length_of_s = (int) floor(L/ds);
123 int half_ls = (int) floor(L/(ds/2));
124
125 double *s = new double[length_of_s];
126 double *sta = new double[length_of_s];
127 double *stb = new double[length_of_s];
128
129 for (i = 0; i<length_of_s; i++) {
130 s[i] = i*ds;
131 sta[i] = 0;
132 stb[i] = 0;
133 }
134
135 /* Convert from s(p) to p(s) and interpolate for accuracy */
136 double *a1x = new double[num_evals];
137 double *a2x = new double[num_evals];
138 double *a3x = new double[num_evals];
139
140 spline(num_evals, 0, 0, 1, 1, s_of_p, sop_num, a1x, a2x, a3x, &iflag);
141
142 double *s_half = new double[half_ls];
143 double *p_of_s_half = new double[half_ls];
144
145 for (i=0; i < half_ls; ++i) {
146 s_half[i] = (double)i*(ds/2);
147 p_of_s_half[i] = seval(num_evals, s_half[i], s_of_p, sop_num, a1x, a2x, a3x, &last);
148 }
149
150 delete[] a1x; delete[] a2x; delete[] a3x;
151 delete[] s_of_p; delete[] sop_num;
152
153 /* Csp is C(s(p)) = [Cx(p(s)) Cy(p(s)) Cz(p(s))] */
154 double *Cspx = new double[length_of_s];
155 double *Cspy = new double[length_of_s];
156 double *Cspz = new double[length_of_s];
157
158 double *p_of_s = new double[length_of_s];
159
160 for (i=0; i<length_of_s; ++i) {
161 p_of_s[i] = p_of_s_half[2*i];
162 Cspx[i] = seval(Lp, p_of_s[i], p, x, c1x, c2x, c3x, &last);
163 Cspy[i] = seval(Lp, p_of_s[i], p, y, c1y, c2y, c3y, &last);
164 Cspz[i] = seval(Lp, p_of_s[i], p, z, c1z, c2z, c3z, &last);
165 }
166 delete[] p_of_s_half;
167
168 /* arrays used by spline function to store coefficients. */
169 double *Csp1x = new double[length_of_s];
170 double *Csp2x = new double[length_of_s];
171 double *Csp3x = new double[length_of_s];
172 double *Csp1y = new double[length_of_s];
173 double *Csp2y = new double[length_of_s];
174 double *Csp3y = new double[length_of_s];
175 double *Csp1z = new double[length_of_s];
176 double *Csp2z = new double[length_of_s];
177 double *Csp3z = new double[length_of_s];
178 spline(length_of_s, 0, 0, 1, 1, s, Cspx, Csp1x, Csp2x, Csp3x, &iflag);
179 spline(length_of_s, 0, 0, 1, 1, s, Cspy, Csp1y, Csp2y, Csp3y, &iflag);
180 spline(length_of_s, 0, 0, 1, 1, s, Cspz, Csp1z, Csp2z, Csp3z, &iflag);
181
182 int size_k = half_ls+2; /* extend of k for RK4 */
183 double *k = new double[size_k]; /* k is the curvature along the curve */
184 for (i=0; i < half_ls; ++i) {
185 double kx = deriv2(length_of_s, s_half[i], s, Csp1x, Csp2x, Csp3x, &last);
186 double ky = deriv2(length_of_s, s_half[i], s, Csp1y, Csp2y, Csp3y, &last);
187 double kz = deriv2(length_of_s, s_half[i], s, Csp1z, Csp2z, Csp3z, &last);
188 k[i] = sqrt(kx*kx + ky*ky + kz*kz); /* the curvature, magnitude of the second derivative of the curve in arc-length parameterization */
189 }
190 k[size_k-2] = k[size_k-3];
191 k[size_k-1] = k[size_k-3];
192
193 delete[] s_half;
194 delete[] Cspx; delete[] Cspy; delete[] Cspz;
195 delete[] Csp1x; delete[] Csp1y; delete[] Csp1z;
196 delete[] Csp2x; delete[] Csp2y; delete[] Csp2z;
197 delete[] Csp3x; delete[] Csp3y; delete[] Csp3z;
198
199 /* computing geomtry dependent constraints (forbidden line curve) */
200 double *sdot = new double[half_ls];
201
202 /* Calculating the upper bound for the time parametrization */
203 /* sdot (which is a non scaled max gradient constaint) as a function of s. */
204 /* sdot is the minimum of gamma*gmax and sqrt(gamma*gmax / k) */
205 double gammagmax = gamma*gmax;
206 for (i=0; i< half_ls; ++i) {
207 double sdot2 = sqrt((gamma*smax) / (fabs(k[i]+(DBL_EPSILON))));
208 if (gammagmax < sdot2)
209 sdot[i] = gammagmax;
210 else
211 sdot[i] = sdot2;
212 }
213
214 double g0gamma = g0*gamma + st0;
215 if (g0gamma < gammagmax)
216 sta[0] = g0gamma;
217 else
218 sta[0] = gammagmax;
219
220 /* Solving ODE Forward */
221
222 for (i=1; i<length_of_s; ++i) {
223 double k_rk[3];
224 k_rk[0] = k[2*i-2];
225 k_rk[1] = k[2*i-1];
226 k_rk[2] = k[2*i];
227
228 double dstds = RungeKutte_riv(ds, sta[i-1], k_rk, smax, gamma);
229 double tmpst = sta[i-1] + dstds;
230
231 if (sdot[2*i+1] < tmpst)
232 sta[i] = sdot[2*i+1];
233 else
234 sta[i] = tmpst;
235 }
236
237 /*Solving ODE Backwards: */
238
239 double max;
240 if(gfin < 0 ) {
241 /*if gfin is not provided */
242 stb[length_of_s-1] = sta[length_of_s - 1];
243 } else {
244
245 if (gfin * gamma > st0)
246 max = gfin*gamma;
247 else
248 max = st0;
249
250 if (gammagmax < max)
251 stb[length_of_s-1] = gammagmax;
252 else
253 stb[length_of_s-1] = max;
254 }
255
256 for (i=length_of_s-2; i>-1; --i) {
257 double k_rk[3];
258 k_rk[0] = k[2*i+2];
259 k_rk[1] = k[2*i+1];
260 k_rk[2] = k[2*i];
261
262 double dstds = RungeKutte_riv(ds, stb[i+1], k_rk, smax, gamma);
263 double tmpst = stb[i+1] + dstds;
264
265 if (sdot[2*i] < tmpst)
266 stb[i] = sdot[2*i];
267 else
268 stb[i] = tmpst;
269 }
270 delete[] k;
271 delete[] sdot;
272
273 /*Final interpolation */
274
275 /* Converting to the time parameterization, t(s) using trapezoidal integration. t(s) = integral (1/st) ds */
276 double st_of_s;
277 if (sta[0] < stb[0])
278 st_of_s = sta[0];
279 else
280 st_of_s = stb[0];
281 double st_ds_i_pre = ds*(1/st_of_s);
282
283 double *t_of_s = new double[length_of_s];
284 t_of_s[0] = 0;
285 for (i=1; i < length_of_s; ++i) {
286 if (sta[i] < stb[i])
287 st_of_s = sta[i];
288 else
289 st_of_s = stb[i];
290
291 double st_ds_i = ds*(1/st_of_s); /* ds * 1/st(s) used in below calculation of t(s) */
292
293 t_of_s[i] = t_of_s[i-1] + (st_ds_i+ st_ds_i_pre)/2;
294
295 st_ds_i_pre = st_ds_i;
296 }
297 delete[] sta; delete[] stb;
298
299 l_t = (int) floor(t_of_s[length_of_s-1]/dt);
300
301
302 /* coefficient arrays for spline interpolation of t(s) to get s(t) */
303 double *t1x = new double[length_of_s];
304 double *t2x = new double[length_of_s];
305 double *t3x = new double[length_of_s];
306
307 spline(length_of_s, 0, 0, 1, 1, t_of_s, s, t1x, t2x, t3x, &iflag);
308
309
310 /* coefficient arrays for spline interpolation of p(s) with s(t) to get p(s(t)) = p(t) */
311 double *p1x = new double[length_of_s];
312 double *p2x = new double[length_of_s];
313 double *p3x = new double[length_of_s];
314
315 spline(length_of_s, 0, 0, 1, 1, s, p_of_s, p1x, p2x, p3x, &iflag);
316
317 double *p_of_t = new double[l_t];
318
319 for (i=0; i < l_t; ++i) {
320 double s_of_t = seval(length_of_s, i*dt, t_of_s, s, t1x, t2x, t3x, &last);
321 p_of_t[i] = seval(length_of_s, s_of_t, s, p_of_s, p1x, p2x, p3x, &last);
322 }
323
324 delete[] t1x; delete[] t2x; delete[] t3x;
325 delete[] t_of_s;
326 delete[] s; delete[] p_of_s;
327 delete[] p1x; delete[] p2x; delete[] p3x;
328
329 /* interpolated k-space trajectory */
330 double *Cx = new double[l_t];
331 double *Cy = new double[l_t];
332 double *Cz = new double[l_t];
333
334 for (i=0; i<l_t; i++) {
335 Cx[i] = seval(Lp, p_of_t[i], p, x, c1x, c2x, c3x, &last);
336 Cy[i] = seval(Lp, p_of_t[i], p, y, c1y, c2y, c3y, &last);
337 Cz[i] = seval(Lp, p_of_t[i], p, z, c1z, c2z, c3z, &last);
338 }
339 delete[] p_of_t; delete[] p;
340 delete[] c1x; delete[] c1y; delete[] c1z;
341 delete[] c2x; delete[] c2y; delete[] c2z;
342 delete[] c3x; delete[] c3y; delete[] c3z;
343
344 /* Final gradient waveforms to be returned */
345 gx = new double[l_t];
346 gy = new double[l_t];
347 gz = new double[l_t];
348
349 for (i=0; i< l_t - 1; ++i) {
350 gx[i] = (Cx[i+1] - Cx[i]) / (gamma * dt);
351 gy[i] = (Cy[i+1] - Cy[i]) / (gamma * dt);
352 gz[i] = (Cz[i+1] - Cz[i]) / (gamma * dt);
353 }
354 delete[] Cx; delete[] Cy; delete[] Cz;
355
356 gx[l_t-1] = 2*gx[l_t-2] - gx[l_t-3];
357 gy[l_t-1] = 2*gy[l_t-2] - gy[l_t-3];
358 gz[l_t-1] = 2*gz[l_t-2] - gz[l_t-3];
359
360}
361
362inline void calcTrajectory(double *gx, double *gy, double *gz, int l_t, double dt, double*& kx, double*& ky, double*& kz, double gamma = 4.25756) {
363 /* k-space trajectory to be returned (calculated by integrating gradient waveforms by trapezoidal integration) */
364 kx = new double[l_t];
365 ky = new double[l_t];
366 kz = new double[l_t];
367
368 double sofarx = 0;
369 double sofary = 0;
370 double sofarz = 0;
371
372 kx[0] = 0.;
373 ky[0] = 0.;
374 kz[0] = 0.;
375
376 for (int i=1; i < l_t; ++i) {
377 sofarx += (gx[i] + gx[i-1]) / 2;
378 sofary += (gy[i] + gy[i-1]) / 2;
379 sofarz += (gz[i] + gz[i-1]) / 2;
380
381 kx[i] = sofarx * dt * gamma;
382 ky[i] = sofary * dt * gamma;
383 kz[i] = sofarz * dt * gamma;
384 }
385}
386
387#endif

-- last change 03.01.2025 | Tony Stoecker | Imprint | Data Protection --