Trajectories for simulated diffusion MRI of turbulent flow in porous media

Dataset contains Lagrangian position and velocity data for a turbulent flow in porous medium. The data consist of 102,170 points over ~0.35 s in increments of 3.5e-4 s. Version 2 Increases the lifetime of the trajectories and greatly improves upon the processing scripts from Version 1. Data are stored in MATLAB .mat format. Scripts are in the MATLAB .m format and were run on R2024b. Functions with names in bold are ones you can run. Non-bolded functions are called by other functions. trajectories*.mat: x,y,z are position data, u,v,w are velocity data. Each file has ~14000 individual trajectories. Separate files for ease of upload/download. trajLoadPos.m/trajLoadVel.m: Scripts that load position or velocity data from the trajectories*.mat files. Note that it will work for any subset of the files, you don't need to download all of them (although it is advisable for good statistics). calcVACF.m: Calculate velocity autocorrelation tensor using trajectory velocity data. Also used to generate a "ground truth" diffusion spectrum tensor, $\boldsymbol{D}(\omega)$ that is derived independently of the simulated MRI measurement. gradientsDOR_v3.m: Generates the double-rotation gradient waveforms used to simulate MRI experiments. See Jiang, et al. Magn. Reson. (2023) for the theory behind this. The gradients are normalised so that the trace of the $b$ tensor is unity. For different $b$ values, multiply the phase distribution by $\sqrt{b}$. simMRI_v4.m simulates multidimensional diffusion MRI experiments on the trajectories. This script produces a "phase tensor", $\boldsymbol{\Phi}$, whose elements are: $\boldsymbol{\Phi} = \int_0^{\tau} \boldsymbol{x}(t)\times\boldsymbol{g}^{\mathrm{T}}(t) \, dt$ This script calclates statistics of the phase tensor elements (variance, kurtosis) to illustrate where the Gaussian phase approximation holds. We save the components of $\boldsymbol{\Phi}$ for each trajectory and each gradient waveform to facilitate further postprocessing in phTensorDists.mat. These are loaded for use by phaseDistTransform.m. phaseDistTransform.m: To obtain a phase distribution $\phi$ realisable by an MRI experiment for any combination of gradient orientation and encoding anisotropy $b_{\mathrm{Mod}}$, we define a transformation tensor, $\boldsymbol{\Gamma}$ that affectuates the reorientation and anisotropy of the gradients. $\boldsymbol{\Gamma} = \begin{bmatrix} R_{xx}b_{\mathrm{Mod},x} & R_{xy}b_{\mathrm{Mod},y} & R_{xz}b_{\mathrm{Mod},z} \\ R_{yx}b_{\mathrm{Mod},x} & R_{yy}b_{\mathrm{Mod},y} & R_{zy}b_{\mathrm{Mod},z} \\ R_{zx}b_{\mathrm{Mod},x} & R_{zy}b_{\mathrm{Mod},y} & R_{zz}b_{\mathrm{Mod},z} \end{bmatrix}\mathrm{,}$ where $R_{ij}$ is the rotation matrix that reorients the gradients from (0,0,1) to the target orientation vector. $b_{\mathrm{Mod}}$ is defined by: $\boldsymbol{b}_{\mathrm{Mod}} = \begin{pmatrix} \sqrt{1-b_{\Delta}} \\ \sqrt{1-b_{\Delta}} \\ \sqrt{1+2b_{\Delta}} \end{pmatrix} $ By taking the tensor product of $\boldsymbol{\Gamma}$ and $\boldsymbol{\Phi}$, we get the phase distribution. This approach means that we can generate a phase distribution without needing to redo the time-consuming time integration step each time. $\phi = \boldsymbol{\Gamma}:\boldsymbol{\Phi}$ Using the phase distribution, we can calculate the NMR signal for a series of $b$ values. Monoexponential decay of $S$ w.r.t. $b$ is expected if the phase distribution is Gaussian and the Stejskal-Tanner model holds. We can then fit the Stejskal-Tanner equation to the NMR signal versus $b$ to find an apparent diffusion coefficient for a given frequency, anisotropy and orientation. These can be compared to the "ground truth" values as a gauge of the method's efficacy. Other MATLAB functions you may need: icosphere.m: Commented out, but useful if you want to generate many different orientation vectors not aligned directly along an axis. https://se.mathworks.com/matlabcentral/fileexchange/50105-icosphere myRotMatFromVecs.m: Turns orientation vector into rotation matrix using the Rodrigues formula. Included here. eigenshuffle.m: Preserves the ordering of the principal components between time and frequency steps. https://se.mathworks.com/matlabcentral/fileexchange/22885-eigenshuffle