OSCR

Laterally Oscillating Trajectory for Undersampling Slices: LOTUS.

Code ↔ Paper

10 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 10 matches · 1 of them tie a paragraph to a whole file, not to given lines: a weak match, whose lines are not tinted
  1. [1] § Theory › G‐Factor Estimation for Non‐Cartesian MRI ↔ operators/mrSampFuncMat.m, the whole file · a weak match · score 0.77 · receiver sensitivity profiles, sampling operator, B0 inhomogeneity, Iterative, masking, MRI
  2. [2] § Methods › Reconstruction ↔ demos/demo_denseMat.m, lines 1–68 · score 0.77 · expanded encoding model, ground truth image, field probe, GB, memory, GPU
  3. [3] § Methods › In Vivo ↔ demos/demo_denseMat.m, lines 1–68 · score 0.73 · field probe system, virtual coils, coil compressed, speed, receivers, noise
  4. [4] § Methods › In Vivo ↔ demos/demo_autoDel.m, lines 1–62 · score 0.71 · field probe system, virtual coils, coil compressed, receivers, noise
  5. [5] § Theory › G‐Factor Estimation for Non‐Cartesian MRI ↔ operators/mrSampFunc.m, lines 1–44 · score 0.65 · receiver sensitivity, domain image, Tikhonov, compensate, Cartesian, operator
  6. [6] § Methods › Reconstruction ↔ computeHarmonics/harmonicsFromRaw.m, lines 1–71 · score 0.60 · encoding model, field probe, coefficients, fitted, iterative, MRI
  7. [7] § Methods › In Vivo ↔ dMRI/nii2kurt.m, lines 1–55 · score 0.55 · matMRI, diffusion tensor, FA, masking
  8. [8] § Methods › Simulations › Investigation of k z Oscillation Period ↔ trajectory/spiralGen.m, lines 90–194 · score 0.55 · Golden ratio, 0–1, oscillation, LOTUS, trajectories
  9. [9] § Theory › Trajectory Design ↔ trajectory/spiralGen.m, lines 1–85 · score 0.54 · gradient relative, Pipe, Nyquist, magnitude, angle, numerical
  10. [10] § Methods › Simulations › Trajectory Comparisons › G‐Factor Maps and Error Metrics ↔ trajectory/spiralGen.m, lines 90–194 · score 0.53 · slew rate, simultaneous slices, msec, resolution, gradient

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

MATLAB · 720 lines · 25 KB · MIT · 3 matches

  1. function [grads,slews,opt] = spiralGen(fovxy,resxy,opt,plotting)
  2. %
  3. % This code is adapted from spiralgen_jgp_12oct.c by Jim Pipe, found at
  4. % https://www.ismrm.org/mri_unbound/sequence.htm
  5. %
  6. % Changes by C. A. Baron:
  7. % - addition of two new spiral types:
  8. % sptype = 4; blipped on z to make multiple spiral planes
  9. % sptype = 5; LOTUS see ISMRM 2025 abstract 1371 (Sothynathan et al)
  10. %
  11. % Inputs:
  12. % fovxy in-plane field of view, in m
  13. % resxy in-plane resolution, in m
  14. % opt options structure. See code for descriptions
  15. %
  17. % Original comments from Pipe:
  18. % /*********************************************
  19. % // Spiral Generation code
  20. % **********************************************
  21. % // Author: Jim Pipe
  22. % // Date: May 2011
  23. % // Rev: Oct 2012
  24. % *********************************************/
  25. % // A Subset of Relevant Literature
  26. % //
  27. % // Spiral Invented:
  28. % // High-speed spiral-scan echo planar NMR imaging-I.
  29. % // Ahn, C.B., Kim, J.H. & Cho, Z.H., IEEE Transactions on Medical Imaging, 5(1) 1986.
  30. % //
  31. % // Spiral Improved:
  32. % // Fast Spiral Coronary Artery Imaging.
  33. % // Meyer CH, Hu BS, Nishimura DG, Macovski A, Magnetic Resonance in Medicine, 28(2) 1992.
  34. % //
  35. % // Variable Density Spiral
  36. % // Reduced aliasing artifacts using variable-density k-space sampling trajectories.
  37. % // Tsai CM, Nishimura DG, Magnetic Resonance in Medicine, 43(3), 2000
  38. % //
  39. % // "SLOPPY" SPIRAL
  40. % // Faster Imaging with Randomly Perturbed Undersampled Spirals and L_1 Reconstruction
  41. % // M. Lustig, J.H. Lee, D.L. Donoho, J.M. Pauly, Proc. of the ISMRM '05
  42. % //
  43. % // FLORET
  44. % // A new design and rationale for 3D orthogonally oversampled k-space trajectories
  45. % // Pipe JG, Zwart NR, Aboussouan EA, Robison RK, Devaraj A, Johnson KO, Mag Res Med 66(5) 2011
  46. % //
  47. % // Distributed Spirals
  48. % // Distributed Spirals: A New Class of 3D k-Space Trajectories
  49. % // Turley D, Pipe JG, Magnetic Resonance in Medicine, in press (also proc of ISMRM '12)
  50. %
  51. % This function
  52. % returns a single spiral arm calculated numerically
  53. %
  54. % The corresponding gradient waveforms are in gxarray and gyarray
  55. % spgrad_na reflects the number of gradient points to reach the end of k-space
  56. % spgrad_nb = spgrad_na + the number of gradient points to ramp G to zero
  57. % spgrad_nc = spgrad_nb + the number of gradient points to rewind k to zero
  58. % spgrad_nd = spgrad_nc + the number of gradient points for first moment compensation
  59. %
  60. % Assignments below indicate units of input parameters
  61. % All units input using kHz, msec, mT, and m!
  62. %
  63. % grad = gm exp(i theta) i.e. gm, theta are magnitude and angle of gradient
  64. % kloc = kr exp(i phi) i.e. kr, phi are magnitude and angle of k-space
  65. % alpha = theta - phi the angle of the gradient relative to that of k-space
  66. % (alpha = Pi/2, you go in a circle
  67. % alpha = 0, you go out radially)
  68. %
  69. % The variable rad_spacing determines the radial spacing
  70. % in units of the Nyquist distance.
  71. % rad_spacing = 1 gives critical sampling
  72. % rad_spacing > 1 gives undersampling
  73. % rad_spacing can vary throughout spiral generation to create variable density spirals
  74. %
  75. % KEY EQUATIONS:
  76. % (1) dkr/dphi = rad_spacing*Nyquist/(2 pi)
  77. % (2) dphi/dt = gamma gm Sin(alpha)/kr
  78. % (3) dkr/dt = gamma gm Cos(alpha)
  79. %
  80. % Solving (1)*(2) = (3) gives
  81. % (4) Tan(alpha) = (2*pi*kr)/(rad_spacing*Nyquist)
  82. %
  83. % *************************************************************/
  84. % /* Initializations */
  85. % /************************************************************/
  87. %%%%% Set limits
  88. maxarray = 128000;
  90. %%%%% Set input defaults
  91. if nargin < 1
  92. fovxy = 0.220/3;
  93. end
  94. if nargin < 2
  95. resxy = 0.002;
  96. end
  97. if nargin < 3
  98. opt = [];
  99. end
  100. if nargin<4 || isempty(plotting)
  101. plotting = true;
  102. end
  103. if ~isfield(opt,'subrast') || isempty(opt.subrast)
  104. opt.subrast = 5; %/* number of numerical cycles per gradient raster time */
  105. end
  106. if ~isfield(opt,'m_dGRast') || isempty(opt.m_dGRast)
  107. opt.m_dGRast = 0.01; % base raster time [msec]
  108. end
  109. if ~isfield(opt,'gamma') || isempty(opt.gamma)
  110. opt.gamma = 42.577; %/* typically 42.577 kHz/mT */
  111. end
  112. if ~isfield(opt,'gmax') || isempty(opt.gmax)
  113. opt.gmax = 30; %/* max gradient amplitude in mT/m */
  114. end
  115. if ~isfield(opt,'slewmax') || isempty(opt.slewmax)
  116. opt.slewmax = 120; %/* max slew rate, in mT/m/msec*/
  117. end
  118. if ~isfield(opt,'gtype') || isempty(opt.gtype)
  119. % 0 = calculate through readout
  120. % 1 = include grad ramp-down
  121. % 2 = include rewinder to end at k=0
  122. % 3 = include first moment comp
  123. opt.gtype = 1;
  124. end
  125. if ~isfield(opt,'fovz') || isempty(opt.fovz)
  126. % Only relevant for opt.sptype == 2
  127. opt.fovz = 0.256; % /* enter in m */
  128. end
  129. if ~isfield(opt,'resz') || isempty(opt.resz)
  130. opt.resz = 0.002; % /* enter in m : this should be true resolution */
  131. end
  132. if ~isfield(opt,'arms') || isempty(opt.arms)
  133. opt.arms = 1; % /* number of spiral interleaves*/
  134. end
  135. if ~isfield(opt,'sptype') || isempty(opt.sptype)
  136. % 0 = Archimedean
  137. % 1 = Cylinder DST
  138. % 2 = Spherical DST
  139. % 3 = Fermat:Floret
  140. % 4 = Archimedean with SMS CAIPI blips (CB 202301)
  141. % 5 = Archimedean with sinusoidal kz for SMS (CB 202301).
  142. % "LOTUS": Laterally Oscillating Trajectory for Undersampling Slices
  143. opt.sptype = 0;
  144. end
  145. % /* the next 4 variables are for variable density spirals */
  146. % /* they create a transition in the radial spacing as the k-space radius goes from 0 to 1, i.e.*/
  147. % /* 0 < kr < us_0 : spacing = Nyquist distance */
  148. % /* us_0 < kr < us_1 : spacing increases to us_r (affected by opt.ustype)*/
  149. % /* us_1 < kr < 1 : spacing = us_r*/
  150. if ~isfield(opt,'ustype') || isempty(opt.ustype)
  151. % rate of change in undersampling
  152. % 0 = linear
  153. % 1 = quadratic
  154. % 2 = hanning */
  155. opt.ustype = 0;
  156. end
  157. if ~isfield(opt,'us_0') || isempty(opt.us_0)
  158. opt.us_0 = 0;
  159. end
  160. if ~isfield(opt,'us_1') || isempty(opt.us_1)
  161. opt.us_1 = 0;
  162. end
  163. if ~isfield(opt,'us_r') || isempty(opt.us_r)
  164. opt.us_r = 1;
  165. end
  166. if ~isfield(opt,'slop_per') || isempty(opt.slop_per)
  167. % For sloppy spirals, this lets us define periodicity in units of iteration loop time */
  168. % set this to zero if you do not want sloppy spirals */
  169. opt.slop_per = 0;
  170. end
  171. % Params for SMS
  172. if ~isfield(opt,'Crate') || isempty(opt.Crate)
  173. opt.Crate = 2; % number of simultaneous slices
  174. end
  175. if ~isfield(opt,'Cdz') || isempty(opt.Cdz)
  176. opt.Cdz = 0.1; % spacing between slices [m]
  177. end
  178. % Params for kz blips for CAIPI-like waveform (sptype == 4)
  179. if ~isfield(opt,'Cf') || isempty(opt.Cf)
  180. opt.Cf = 0.3; % fractional slew rate reserved for blips. 0.2 or 0.3 seems to be a good tradeoff
  181. end
  182. if ~isfield(opt,'CstaggerBlip') || isempty(opt.CstaggerBlip)
  183. opt.CstaggerBlip = 0; %// to stagger blips to try to have fewer k-space gaps
  184. end
  185. % Params for kz oscillations for CAIPI-like waveform (sptype == 5)
  186. if ~isfield(opt,'CphiFact') || isempty(opt.CphiFact)
  187. % period of kz wave is 1/CphiFact larger than kxy rotation period.
  188. % Good choices take many repetitions to get back to an integer to
  189. % moreorless isotropically sample kz. It may also be a good idea to
  190. % choose a value close to 1 so that all channels have similar freq
  191. % content, which makes it easier to avoid vibrational
  192. % resonances.
  193. opt.CphiFact = 0.618; % Equal to 1/golden ratio.
  194. end
  197. %%%%% Internal variable calculation
  198. rast = opt.m_dGRast / opt.subrast; %/* calculation "raster time" in msec */
  199. if ( (opt.sptype >= 4) && (opt.Crate < 2) )
  200. % If only one slice, no need to account for SMS
  201. opt.sptype = 0;
  202. end
  204. %%%%% Error checking
  205. if (opt.CstaggerBlip)
  206. error('Proper handling of rewinding time not tested for staggered blips')
  207. end
  209. %%%%% Start computations
  210. nyquist = opt.arms/fovxy; %/* radial distance per arm to meet the Nyquist limit*/
  211. gamrast = opt.gamma*rast; %/* gamrast*g = dk*/
  212. dgc = opt.slewmax*rast; %/* the most the gradients can change in 1 raster period*/
  213. sub_gamrast = opt.subrast*gamrast;
  214. sub_dgc = opt.subrast*dgc;
  217. uz=0;
  218. gx=0;
  219. gy=0;
  220. gz=0;
  221. kx = zeros(opt.subrast*maxarray,1);
  222. ky = zeros(opt.subrast*maxarray,1);
  223. kz = zeros(opt.subrast*maxarray,1);
  224. gsign = ones(opt.subrast*maxarray,1);
  225. gxarray = zeros(maxarray,1);
  226. gyarray = zeros(maxarray,1);
  227. gzarray = zeros(maxarray,1);
  229. krmax = 0.5/resxy;
  230. kzmax = 0.5/opt.resz;
  231. krmax2 = krmax*krmax;
  232. kzmax2 = kzmax*kzmax;
  233. krlim = krmax*(1.-(resxy/fovxy));
  235. if (opt.sptype==4)
  236. %/* Determine k-space spacing based on "FOV" Crate*dz. This is the
  237. % typical dkz provided by a blip.
  238. CkbTot = 1.0 / (opt.Crate*opt.Cdz);
  239. %/* Determine total k-space step provided by largest blip, which
  240. % brings you back to the starting k-space position. We design to
  241. % this, then scale the smaller blips back down. */
  242. CkbTot = (opt.Crate - 1.0) * CkbTot;
  243. %// Determine duration based on basic k and grad relations [ms]
  244. Ct = 2*sqrt(CkbTot / (opt.gamma * opt.Cf * opt.slewmax));
  245. %// Find total number of points. Make total time a multiple of the twice the raster time
  246. Cnp = floor(Ct/(2.0*opt.m_dGRast) + 1.0) * 2 * opt.subrast;
  247. Ct = Cnp * rast;
  248. %// Find gradient change per point in array using k to grad relationship
  249. Cgmax = 2.0 * CkbTot / (opt.gamma * Ct);
  250. CGstep = Cgmax / (Cnp/2);
  251. %// Fill array with k-space values
  252. Cgb = zeros(Cnp,1);
  253. Cgb(1) = 0;
  254. for i=2:Cnp %(i=1;i<Cnp;i++)
  255. if (i<=Cnp/2)
  256. Cgb(i) = (i-1)*CGstep;
  257. else
  258. Cgb(i) = (Cnp-i+1)*CGstep;
  259. end
  260. end
  261. Cbstart = -1;
  262. CbstartPrev = -1;
  263. Ccurrblip = -1;
  264. CcurrblipPrev = -1;
  265. Cadd = 0; %// To stagger where the blips start a bit
  266. elseif (opt.sptype == 5)
  267. %// Determine k-space spacing based on "FOV" Crate*dz;
  268. CkbTot = 1.0 / (opt.Crate*opt.Cdz);
  269. %// Determine total span to go from min to max k
  270. CkbTot = (opt.Crate - 1.0) * CkbTot;
  271. isStart = 1;
  272. isStartRecord = ones(opt.subrast*maxarray,1);
  273. phiUnwrapped = zeros(opt.subrast*maxarray,1);
  274. end
  276. %/* start out spiral going radially at max slew-rate for 2 time-points */
  277. kx(1) = 0;
  278. ky(1) = 0;
  279. kx(2) = gamrast*dgc;
  280. ky(2) = 0;
  281. kx(3) = 3*gamrast*dgc;
  282. ky(3) = 0;
  284. %// IF SPHERE
  285. if (opt.sptype == 2)
  286. kz(1) = kzmax;
  287. kz(2) = sqrt(kzmax2*(1-((kx(1)*kx(1)+ky(1)*ky(1))/krmax2))); %// stay on surface of ellipsoid
  288. kz(3) = sqrt(kzmax2*(1-((kx(2)*kx(2)+ky(2)*ky(2))/krmax2))); %// stay on surface of ellipsoid
  289. end
  291. nLoops = 0;
  292. i = 3;
  293. kr = kx(3);
  295. % /******************************/
  296. % /* LOOP UNTIL YOU HIT MAX RES */
  297. % /******************************/
  298. while ((kr <= krlim) && (i < opt.subrast*maxarray-1) )
  299. if (nLoops > 10*opt.subrast*maxarray)
  300. error('Spiral gen failure')
  301. end
  302. if (i<2)
  303. error('Spiral gen failure on rewinding time')
  304. end
  306. % /**************************/
  307. % /*** STEP 1: Determine the direction (ux,uy) of the gradient at ~(i+0.5) */
  308. % /**************************/
  309. % /* calculate dk/rast = opt.gamma G*/
  310. kmx = 1.5*kx(i) - 0.5*kx(i-1);
  311. kmy = 1.5*ky(i) - 0.5*ky(i-1);
  312. kmr = sqrt(kmx*kmx + kmy*kmy);
  314. % /////////////////////////////
  315. % // Start rad_spacing logic //
  316. % /////////////////////////////
  317. rnorm = 2*resxy*kmr; %/* the k-space radius, normalized to go from 0 to 1 */
  319. %/* determine the undersample factor */
  320. if (rnorm <= opt.us_0)
  321. rad_spacing = 1;
  322. elseif (rnorm < opt.us_1)
  323. us_i = (rnorm-opt.us_0)/(opt.us_1 - opt.us_0); %/* goes from 0 to 1 as rnorm goes from us_0 to us_1*/
  324. if (opt.ustype == 0)
  325. %/* linearly changing undersampling*/
  326. rad_spacing = 1. + (opt.us_r - 1.)*us_i;
  327. elseif (opt.ustype == 1)
  328. %/* quadratically changing undersampling*/
  329. rad_spacing = 1. + (opt.us_r - 1.)*us_i*us_i;
  330. elseif (opt.ustype == 2)
  331. %/* Hanning-type change in undersampling */
  332. rad_spacing = 1. + (opt.us_r - 1.)*0.5*(1.-cos(us_i*M_PI));
  333. end
  334. else
  335. rad_spacing = opt.us_r;
  336. end
  338. %/* Undersample spiral for Spherical-Distributed Spiral */
  339. if (opt.sptype == 2)
  340. if (rnorm < 1.0)
  341. rad_spacing = min(opt.fovz/opt.resz, rad_spacing/sqrt(1.0 - (rnorm*rnorm)));
  342. else
  343. rad_spacing = opt.fovz/opt.resz;
  344. end
  345. end
  347. %/* MAKE FERMAT SPIRAL FOR FLORET*/
  348. if (opt.sptype == 3 && rnorm > 0)
  349. rad_spacing = rad_spacing / rnorm;
  350. end
  352. %/* Sloppy Spirals - add variability to rad_spacing for reduced aliasing coherence */
  353. % // A couple different options here are commented out
  354. % // Lots of ways to be sloppy
  355. if (opt.slop_per > 0)
  356. % // rad_spacing = MAX(1., (rad_spacing + ((rad_spacing-1.)*sin(2.*M_PI*(double)(i)/opt.slop_per))));
  357. % // rad_spacing += (rad_spacing-1.)*sin(2.*M_PI*opt.slop_per*atan2(ky(i),kx(i)));
  358. rad_spacing = rad_spacing + (rad_spacing-1)*sin(2*pi*opt.slop_per*rnorm);
  359. end
  361. % ///////////////////////////
  362. % // End rad_spacing logic //
  363. % ///////////////////////////
  365. %/* See the Key Equation 4 at the beginning of the code */
  366. alpha = atan(2*pi*kmr/(rad_spacing*nyquist));
  367. phi = atan2(kmy,kmx);
  368. theta = phi + alpha;
  370. ux = cos(theta);
  371. uy = sin(theta);
  373. % // IF SPHERICAL DST
  374. % // u dot km is zero if moving on a sphere (km is radial, u is tangential,
  375. % // thus km stays on the sphere)
  376. % // We are on an ellipsoid, but can normalize u and km by krmax and kzmax to make this work
  377. % // The final gradient vector (ux uy uz) will be tangential to the sphere
  378. if (opt.sptype == 2)
  379. kmz = 1.5*kz(i) - 0.5*kz(i-1);
  380. uz = -((ux*kmx + uy*kmy)/krmax2)*(kzmax2/kmz);
  381. umag = sqrt(ux*ux + uy*uy + uz*uz);
  382. ux = ux/umag;
  383. uy = uy/umag;
  384. uz = uz/umag;
  385. gz = (kz(i) - kz(i-1))/gamrast;
  386. elseif (opt.sptype == 4)
  387. % Set gradient for pre-defined CAIPI blip
  388. gz = (kz(i) - kz(i-1))/gamrast;
  389. %// Start a blip whenever kx goes from pos to neg
  390. if ( (kx(i-1) > 0) && (kx(i) < 0) )
  391. if ( (Cbstart > 0) && (i >= Cbstart) && (i < Cbstart + Cnp ) )
  392. error('CAIPI blips overlapping')
  393. end
  394. CbstartPrev = Cbstart;
  395. Cbstart = i + Cadd;
  396. CcurrblipPrev = Ccurrblip;
  397. Ccurrblip = Ccurrblip + 1;
  398. if (Ccurrblip > opt.Crate)
  399. Ccurrblip = 1;
  400. end
  401. if opt.CstaggerBlip
  402. Cadd = Cadd + round(Cnp/opt.Crate);
  403. if Cadd > Cnp
  404. Cadd = 0;
  405. end
  406. end
  407. end
  408. %// Check if in blip
  409. if ( (Cbstart > 0) && (i >= Cbstart) && (i < Cbstart + Cnp ) )
  410. %// Scale blip and set polarity
  411. if (Ccurrblip == 0)
  412. %// First blip is scaled differently to make kz symmetric
  413. Cfact = -0.5;
  414. elseif (Ccurrblip < opt.Crate)
  415. Cfact = 1/(opt.Crate-1);
  416. else
  417. Cfact = -1;
  418. end
  419. gznext = Cfact*Cgb(i-Cbstart+1);
  420. else
  421. gznext = 0.0;
  422. end
  423. elseif (opt.sptype == 5)
  424. %/* Find unwrapped dphi. */
  425. dphi = atan2(ky(i),kx(i)) - atan2(ky(i-1),kx(i-1));
  426. if abs(dphi)>3*pi/4
  427. if (dphi<0)
  428. dphi = dphi + 2*pi;
  429. else
  430. dphi = dphi - 2*pi;
  431. end
  432. end
  433. phiUnwrapped(i) = phiUnwrapped(i-1) + dphi;
  434. % Our target is kz = CkbTot/2*cos(opt.CphiFact*phi).
  435. % Thus, dk/dphi = CkbTot/2*opt.CphiFact*sin(opt.CphiFact*phi)
  436. % From comments at top of file, dphi/dt = gamma*Gxy*sin(alpha)/kmr
  437. % Multiplying these equations: dk/dt = gamma*Gxy*sin(alpha)/kmr*CkbTot/2*opt.CphiFact*sin(opt.CphiFact*phi)
  438. % Recognizing that gamma*Gz = dk/dt, using ux,uy,uz for Gx,Gy,Gz, and rearranging yields:
  439. uz = sqrt(ux^2+uy^2)*opt.CphiFact*CkbTot/2/kmr*...
  440. sin(opt.CphiFact*phiUnwrapped(i))*sin(alpha);
  441. % If we just use the above uz, the kz span will not be centered.
  442. % So, we use the first half-period of cos(opt.CphiFact*phi) as a
  443. % "prephasor" to get to the edge of the desired span. We can do
  444. % this by scaling uz by 0.5 during the first pi radians of opt.CphiFact*phi
  445. isStartRecord(i) = isStart;
  446. if (isStart)
  447. uz = 0.5*uz;
  448. if (opt.CphiFact*phiUnwrapped(i) >= pi)
  449. isStart = 0;
  450. end
  451. end
  452. %// Normalize the unit vector
  453. umag = sqrt(ux*ux + uy*uy + uz*uz);
  454. ux = ux/umag;
  455. uy = uy/umag;
  456. uz = uz/umag;
  457. gz = (kz(i) - kz(i-1))/gamrast;
  458. end
  460. % /**************************/
  461. % /*** STEP 2: Find largest gradient magnitude with available slew */
  462. % /**************************/
  464. %/* Current gradient*/
  465. gx = (kx(i) - kx(i-1))/gamrast;
  466. gy = (ky(i) - ky(i-1))/gamrast;
  468. % /*
  469. % // solve for gm using the quadratic equation |gm u - g| = dgc
  470. % // which is
  471. % // (gm u - g)(gm u* - g*) = dgc^2
  472. % // which gives
  473. % // gm^2 (u u*) - gm (g u* + u g*) + g g* - dgc^2 = 0
  474. %
  475. % // Replacing u u* with 1 (i.e. u is a unit vector) and
  476. % // replacing (g u* + u g*) with 2 Real[g u*]
  477. % // this is
  478. % // gm^2 + gm (2 b) + c = 0
  479. % // giving
  480. % // gm = -b +/- Sqrt(b^2 - c)
  481. % // The variable "term" = (b^2 - c) will be positive if we can meet the desired new gradient
  482. % */
  483. if (opt.sptype == 4)
  484. %// Only allow slew not reserved for blips. Ignore gz
  485. %// keep slew reduced for whole waveform. Could be more efficient, but this is easy
  486. term = dgc*dgc*(1-opt.Cf^2) - (gx*gx + gy*gy) + (ux*gx + uy*gy)*(ux*gx + uy*gy);
  487. else
  488. term = dgc*dgc - (gx*gx + gy*gy + gz*gz) + (ux*gx + uy*gy + uz*gz)*(ux*gx + uy*gy + uz*gz);
  489. end
  491. if (term >= 0)
  492. % // Slew constraint is met! Now assign next gradient and then next k value
  493. % // NOTE gsign is +1 or -1
  494. % // if gsign is positive, we are using slew to speed up (increase gm) as much as possible
  495. % // if gsign is negative, we are using slew to slow down (decrease gm) as much as possible
  496. if (opt.sptype == 4)
  497. %// Account for fixed blip gradient contributing to net grad
  498. %// keep xy max grad reduced for whole waveform. Could be more efficient, but this is easy
  499. gm = min((ux*gx + uy*gy) + gsign(i)*sqrt(term),sqrt(opt.gmax^2-Cgmax^2));
  500. else
  501. gm = min((ux*gx + uy*gy + uz*gz) + gsign(i)*sqrt(term),opt.gmax);
  502. end
  503. gx = gm*ux;
  504. gy = gm*uy;
  506. kx(i+1) = kx(i) + gx*gamrast;
  507. ky(i+1) = ky(i) + gy*gamrast;
  509. %// If SPHERE
  510. if (opt.sptype == 2)
  511. kz(i+1) = sqrt(kzmax2*(1.-((kx(i+1)*kx(i+1)+ky(i+1)*ky(i+1))/krmax2))); %// stay on surface of ellipsoid
  512. elseif (opt.sptype == 4)
  513. kz(i+1) = kz(i) + gznext*gamrast;
  514. elseif (opt.sptype == 5)
  515. gz = gm*uz;
  516. kz(i+1) = kz(i) + gz*gamrast;
  517. end
  519. i = i+1;
  520. else
  521. % // We can't go further without violating the slew rate
  522. % // This means that we've sped up too fast to turn here at the desired curvature
  523. % // We are going to iteratively go back in time and slow down, rather than speed up, at max slew
  524. % // Here we'll keep looking back until gsign is positive, then add another negative gsign, just far enough to make the current corner
  525. while ((i>4) && (gsign(i-1) == -1))
  526. i = i-1;
  527. end
  528. gsign(i-1) = -1;
  529. i = i-2;
  530. if (opt.sptype == 4) && (i<=Cbstart)
  531. % We rewound past a blip start
  532. Cbstart = CbstartPrev;
  533. Ccurrblip = CcurrblipPrev;
  534. elseif (opt.sptype == 5)
  535. isStart = isStartRecord(i);
  536. end
  537. end
  539. kr = sqrt(kx(i)*kx(i) + ky(i)*ky(i));
  541. nLoops = nLoops + 1;
  543. end % End main kr while loop
  545. i_end = i;
  547. % //********************************************
  548. % // DONE LOOPING FOR SAMPLING PORTION
  549. % // recast k to g while subsampling by opt.subrast
  550. % //********************************************
  551. % TODO: vectorize this
  552. gxsum = 0;
  553. gysum = 0;
  554. gzsum = 0;
  555. for j = 1:floor(i_end/opt.subrast) %(j=1;j<=(i_end/opt.subrast);j++)
  556. i1 = j*opt.subrast + 1;
  557. i0 = (j-1)*opt.subrast + 1;
  558. gxarray(j) = ( (kx(i1)-kx(i0))/sub_gamrast );
  559. gyarray(j) = ( (ky(i1)-ky(i0))/sub_gamrast );
  560. gzarray(j) = ( (kz(i1)-kz(i0))/sub_gamrast );
  561. gxsum = gxsum + gxarray(j);
  562. gysum = gysum + gyarray(j);
  563. gzsum = gzsum + gzarray(j);
  564. end
  565. spgrad_na = j;
  567. %// recalculate these ending gradient points
  568. gm = sqrt(gxarray(spgrad_na-1)*gxarray(spgrad_na-1) +...
  569. gyarray(spgrad_na-1)*gyarray(spgrad_na-1) +...
  570. gzarray(spgrad_na-1)*gzarray(spgrad_na-1));
  571. ux = gxarray(spgrad_na-1)/gm;
  572. uy = gyarray(spgrad_na-1)/gm;
  573. uz = gzarray(spgrad_na-1)/gm;
  575. % //**************************************************
  576. % // NOW, if requested via gtype, go to g=0 and k=0
  577. % // I've tried other ways to be faster, can't find them
  578. % //**************************************************
  580. % // first we'll ramp gradients to zero
  581. % // note {ux,uy} is still pointing in the gradient direction
  582. % TODO: vectorize
  583. if (opt.gtype > 0)
  584. gz_sum_ramp = 0;
  585. while ((gm > 0) && (j < maxarray))
  586. gm = max(0,gm - sub_dgc);
  587. gxarray(j) = gm*ux;
  588. gyarray(j) = gm*uy;
  589. gzarray(j) = gm*uz;
  590. gxsum = gxsum + gxarray(j);
  591. gysum = gysum + gyarray(j);
  592. gzsum = gzsum + gzarray(j);
  593. gz_sum_ramp = gz_sum_ramp + gzarray(j);
  594. j = j+1;
  595. end
  596. end
  597. spgrad_nb = j;
  599. % // now point gradient towards the k-space origin
  600. % // {ux,uy} will be a unit vector in that direction
  601. if (opt.gtype > 1)
  602. % /* NOTE: spherical needs a prephaser not a rewinder
  603. % * so just rewind x and y in that case */
  604. gsum = sqrt(gxsum*gxsum + gysum*gysum + gzsum*gzsum);
  605. if (opt.sptype == 2 )
  606. gsum = sqrt(gxsum*gxsum + gysum*gysum + gz_sum_ramp*gz_sum_ramp);
  607. end
  608. gsum0 = gsum;
  609. ux = -gxsum/gsum;
  610. uy = -gysum/gsum;
  611. uz = -gzsum/gsum;
  612. if (opt.sptype == 2)
  613. uz = -gz_sum_ramp/gsum;
  614. end
  615. gsum_ramp = 0.5*gm*(gm/sub_dgc); %/* this is *roughly* how much the area changes if we ramp down the gradient NOW*/
  616. %/* this value is zero right now (gm = 0), but it will make sense below */
  617. %// increase gm while we can
  618. while ((gsum_ramp < gsum) && (j < maxarray))
  619. gm = min(opt.gmax,gm+sub_dgc);
  620. gxarray(j) = gm*ux;
  621. gyarray(j) = gm*uy;
  622. gzarray(j) = gm*uz;
  623. gsum = gsum - gm;
  624. j = j+1;
  625. gsum_ramp = 0.5*gm*(gm/sub_dgc); %/* see - now this makes sense; this tells us when to start ramping down */
  626. end
  628. % // We've overshot it by a tiny bit, but we'll fix that later
  629. % // Ramp down for now
  630. while ((gm > 0) && (j < maxarray))
  631. gm = max(0,gm-sub_dgc);
  632. gxarray(j) = gm*ux;
  633. gyarray(j) = gm*uy;
  634. gzarray(j) = gm*uz;
  635. gsum = gsum - gm;
  636. j = j+1;
  637. end
  638. spgrad_nc = j;
  640. %// OK - gm is zero, but gsum is probably not EXACTLY zero. Now scale the rewinder to make the sum exactly zero
  641. gradtweak = gsum0/(gsum0-gsum);
  642. for j = spgrad_nb+1:spgrad_nc %(j=(spgrad_nb); j<(spgrad_nc); j++)
  643. gxarray(j) = (gradtweak)*gxarray(j);
  644. gyarray(j) = (gradtweak)*gyarray(j);
  645. gzarray(j) = (gradtweak)*gzarray(j);
  646. end
  647. end
  649. gxarray = gxarray(1:j);
  650. gyarray = gyarray(1:j);
  651. gzarray = gzarray(1:j);
  653. grads = cat(2, gxarray, gyarray, gzarray);
  654. slews = diff(grads,1,1)/opt.m_dGRast;
  657. %% Plotting
  658. if plotting
  659. figure;
  660. n1 = 2;
  661. n2 = 4;
  662. subplot(n1,n2,1);
  663. plot(grads);
  664. title('gradients')
  666. subplot(n1,n2,2)
  667. kx = cumsum(gxarray)*opt.gamma*opt.m_dGRast;
  668. ky = cumsum(gyarray)*opt.gamma*opt.m_dGRast;
  669. kz = cumsum(gzarray*opt.gamma*opt.m_dGRast);
  670. plot(kx, ky)
  671. if opt.sptype == 5
  672. zc = [abs(diff(sign(kz)))>0;false];
  673. kx_a = kx(zc);
  674. ky_a = ky(zc);
  675. hold('all');
  676. plot(kx_a,ky_a,'o')
  677. end
  678. title('k-space traj xy with kz zero crossings')
  680. subplot(n1,n2,3)
  681. plot3(cumsum(gxarray*opt.gamma*opt.m_dGRast), cumsum(gyarray*opt.gamma*opt.m_dGRast), cumsum(gzarray*opt.gamma*opt.m_dGRast))
  682. title('k-space traj xyz')
  684. subplot(n1,n2,4)
  685. % plot(cumsum(gxarray*opt.gamma*opt.m_dGRast))
  686. % hold('all')
  687. % plot(cumsum(gyarray*opt.gamma*opt.m_dGRast))
  688. plot(cumsum(gzarray*opt.gamma*opt.m_dGRast))
  689. title('k-space traj z')
  691. subplot(n1,n2,n2+1);
  692. plot(sqrt(sum(slews.^2,2)));
  693. title('net slew')
  695. subplot(n1,n2,n2+2);
  696. plot(sqrt(sum(grads.^2,2)));
  697. title('net grad')
  699. % Freq analysis notes:
  700. % Gmax and slew has the largest effect.
  701. % Undersampling and variable density has little effect
  702. % CAIPI grads are negligible compared to others.
  703. subplot(n1,n2,n2+3)
  704. fmaxp = 3000;
  705. Nf = 10*length(gxarray);
  706. fmax = 1000*0.5/opt.m_dGRast; % Hz
  707. psd = fft([gxarray,gyarray,gzarray],Nf,1);
  708. psd = psd(1:ceil(Nf/2),:).*conj(psd(1:ceil(Nf/2),:));
  709. f = linspace(0,fmax,size(psd,1))';
  710. plot(f,psd)
  711. xlim([0,fmaxp])
  712. cent = sum(f.*psd,1)./sum(psd,1);
  713. [~, pk] = max(psd(:,1)); pk = f(pk);
  714. bw = find(psd(:,1) > max(psd(:,1))/2);
  715. bw = f(bw(end)) - f(bw(1));
  716. text(0.5*fmaxp, 0.9*max(psd(:,1)),...
  717. sprintf('cent = %d Hz\npeak = %d Hz\nbw = %d Hz\nrate %d\nvd %.2f',round(cent(1)),round(pk(1)),round(bw),opt.us_r,opt.us_1))
  718. end

spiralGen.m at commit 17825cc, under MIT · at the source

Overview

Authors: Mayuri Sothynathan1,2, Paul I. Dubovan3,4, Corey A. Baron1,5
ORCID iDs: Mayuri Sothynathan, Corey A. Baron
  1. Centre for Functional and Metabolic Mapping (CFMM), Robarts Research Institute, Western University London Ontario Canada
  2. Department of Biomedical Engineering Faculty of Engineering, Western University London Ontario Canada
  3. Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital Charlestown Massachusetts USA
  4. Department of Radiology Harvard Medical School Boston Massachusetts USA
  5. Department of Medical Biophysics Schulich School of Medicine and Dentistry, Western University London Ontario Canada
Institutions: Western University (Canada); Robarts Research Institute (Canada); Harvard University (United States); Massachusetts General Hospital (United States); Athinoula A. Martinos Center for Biomedical Imaging (United States)
Journal: Magnetic resonance in medicine, volume 96, issue 4, pages 1682-1695
Dates: received 4 February 2026; accepted 1 June 2026; published online 9 June 2026; in print October 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1002/mrm.70469 · PMID 42265901 · PMCID PMC13421052 · OpenAlex W7164123668
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism), computational (subfield)
Methods: Machine learning, fMRI & imaging
Keywords: diffusion MRI, g‐factor, non‐Cartesian, simultaneous multislice, spiral
MeSH: Brain*, Diffusion Magnetic Resonance Imaging*, Image Processing, Computer-Assisted*, Imaging, Three-Dimensional*, Algorithms, Anisotropy, Computer Simulation, Humans, Image Enhancement, Image Interpretation, Computer-Assisted, Phantoms, Imaging, Reproducibility of Results, Signal-To-Noise Ratio (* major topic)
Journal subjects: Imaging Methodology
Topic: Advanced Neuroimaging Techniques and Applications (Radiology, Nuclear Medicine and Imaging, Medicine), according to OpenAlex
Funding: Ontario Graduate Scholarship; NSERC Discovery Grants (RGPIN‐2025‐05597)
Citations: not cited yet (Europe PMC); 42 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

Its files are read in the Code ↔ Paper reader above, with 10 matches between paragraphs and lines of code.

cfmm/matlab/matmri

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 17825cc25a717696e95e921018e7487d1c156fbf, 10 August 2026
Languages: MATLAB (65)
Size: 76 files, 65 scripts
Software Heritage: not archived
Found in: “Data Availability Statement”
Holds: README, license file, tests
Not found: CITATION.cff, environment file, continuous integration, documentation
Tools: Parallel Computing Toolbox (12 files), Image Processing Toolbox (6 files), Optimization Toolbox (1 file), Statistics and Machine Learning Toolbox (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
67 files
At the source: gitlab.com/cfmm/matlab/matmri

Zenodo 4495476

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
At the source:

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 65 scripts, each with its path and the digest of its content;
  • 10 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1002/mrm.70469.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

  • Publisher: n/a → Wiley

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 3 authors, 5 keywords, 13 MeSH terms, 2 funders, 38 references.

Cite

This paper

Sothynathan, M., Dubovan, P. I., & Baron, C. A. (2026). Laterally Oscillating Trajectory for Undersampling Slices: LOTUS. Magnetic resonance in medicine, 96(4), 1682-1695. https://doi.org/10.1002/mrm.70469

BibTeX

@article{sothynathan2026laterally,
author = {Sothynathan, Mayuri and Dubovan, Paul I. and Baron, Corey A.},
title = {{Laterally Oscillating Trajectory for Undersampling Slices: LOTUS}},
journal = {Magnetic resonance in medicine},
year = {2026},
month = jun,
volume = {96},
number = {4},
pages = {1682--1695},
publisher = {Wiley},
issn = {0740-3194},
doi = {10.1002/mrm.70469},
url = {https://doi.org/10.1002/mrm.70469},
pmid = {42265901},
pmcid = {PMC13421052}
}

RIS

TY - JOUR
AU - Sothynathan, Mayuri
AU - Dubovan, Paul I.
AU - Baron, Corey A.
TI - Laterally Oscillating Trajectory for Undersampling Slices: LOTUS
T2 - Magnetic resonance in medicine
J2 - Magn Reson Med
PY - 2026
DA - 2026/06/09
VL - 96
IS - 4
SP - 1682
EP - 1695
SN - 0740-3194
PB - Wiley
DO - 10.1002/mrm.70469
UR - https://doi.org/10.1002/mrm.70469
LA - en
ER -

CSL-JSON

{
"id": "10.1002/mrm.70469",
"type": "article-journal",
"title": "Laterally Oscillating Trajectory for Undersampling Slices: LOTUS",
"container-title": "Magnetic resonance in medicine",
"author": [
{
"family": "Sothynathan",
"given": "Mayuri"
},
{
"family": "Dubovan",
"given": "Paul I."
},
{
"family": "Baron",
"given": "Corey A."
}
],
"container-title-short": "Magn Reson Med",
"volume": "96",
"issue": "4",
"page": "1682-1695",
"DOI": "10.1002/mrm.70469",
"PMID": "42265901",
"PMCID": "PMC13421052",
"ISSN": "0740-3194",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/mrm.70469",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
9
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1002/mrm.70439 [code]
Group-Patch Joint Compression: Compressing Dynamic B&lt;sub&gt;0&lt;/sub&gt; and Static RF Spatial Modulations Across k-Space Subregion Groups for Highly Accelerated MRI.
Journal: Magnetic resonance in medicine
In common: Image Processing Toolbox, structural MRI / diffusion, 8 references
[2] doi:10.1162/imag.a.1262 [code]
Frame-wise multi-echo distortion correction for superior functional MRI.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Optimization Toolbox, Parallel Computing Toolbox, Image Processing Toolbox, 1 other tool, 1 reference
[3] doi:10.1038/s41467-026-74215-5 [code]
Multi-metric evaluations of acute psychedelic effects on fMRI brain entropy.
Journal: Nature communications
In common: Optimization Toolbox, Parallel Computing Toolbox, Image Processing Toolbox, 1 other tool, computational
[4] doi:10.1002/mrm.70433 [code]
Single-Shot 2D Radial Echo Planar Imaging for Functional MRI.
Journal: Magnetic resonance in medicine
In common: 4 references
[5] doi:10.1002/hbm.70553 [code]
Axon Diameter Mapping in the Living Human Brain with Ultra-High-Gradient Diffusion MRI at 500 mT/m Gradient Strength.
Journal: Human brain mapping
In common: Optimization Toolbox, Parallel Computing Toolbox, Image Processing Toolbox, 1 other tool, structural MRI / diffusion
[6] doi:10.1038/s42003-026-10276-y [code]
The cellular correlates and adolescent reorganisation of cortical myelination networks in the common marmoset.
Journal: Communications biology
In common: Optimization Toolbox, Parallel Computing Toolbox, Image Processing Toolbox, 1 other tool, structural MRI / diffusion
[7] doi:10.1002/mrm.70366 [code]
Mesoscale Whole-Brain T&lt;sub&gt;2&lt;/sub&gt;*-Weighted and Associated Quantitative MRI in Humans at 10.5 T.
Journal: Magnetic resonance in medicine
In common: Image Processing Toolbox, Statistics and Machine Learning Toolbox, structural MRI / diffusion, 2 references
[8] doi:10.1158/2767-9764.crc-25-0579 [code]
Glioblastoma Subtypes Exhibit Distinct Migration Mechanics and Immune Responses.
Journal: Cancer research communications
In common: Optimization Toolbox, Parallel Computing Toolbox, Image Processing Toolbox, 1 other tool
[9] doi:10.1371/journal.pcbi.1014563 [code]
Single pulse electrical stimulation in white matter modulates iEEG visual responses in human early visual cortex.
Journal: PLoS computational biology
In common: Optimization Toolbox, Parallel Computing Toolbox, Image Processing Toolbox, 1 other tool
[10] doi:10.1007/s00429-026-03152-2 [code]
Limb-selective regions in the lateral temporal lobe shrink from childhood to adulthood.
Journal: Brain structure & function
In common: Optimization Toolbox, Parallel Computing Toolbox, Image Processing Toolbox, 1 other tool

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.